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ABSTRACT 


For the optimal use of high precision Lunar Laser Ranging (LLR), an , 
investigation regarding a clear definition of the underlying coordinate systems, 
identification of estimable quantities, favorable station geometry and optimal 
observation schedule is given. 

hi Section 2, the least squares adjustment formulation for range- 
differencing is presented. Taking advantage of the earth -moon geometry, this 
procedure determines the coordinate differences of the stations particularly 
well. The body-fixed motions of the celestial pole (polar motion) and the earth 
rotation parameter are derived from an orthogonal transformation relative to a 
standard epoch. This is accomplished by a second least squares solution which 
utilizes the estimable parameters of the first adjustment as new "observations. 
A separation between earth rotation variations and ephemeris errors in lunar 
right ascension is not possible. Various station distributions are analyzed. A 
station geometry consisting of two north-south lines, being separated in longi- 
tude by 90°, and one east-west line determine the three orientation parameters 
virtually independent of ephemeris errors in declination. However, they 
include the common motions of the stations due to crustal motions . ' 

The third section presents various analyses of variance models and 
numerical results . The simplifications consist of neglecting the earth rotation 
during the travel time of the pulse, hi some models, the terms of the charac- 
teristic order of 1/60, i.e. , those terms depending on the ratio of geocentric 
station distance to geocentric reflector distance, are neglected and the declina- 
tion is taken constant during one interval. The analysis shows that for the 
given station distribution and an adequate observation schedule, the orientation 
parameters can be given daily with at least the measurement accuracy. 
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LIST OF SYMBOLS 


This list of symbols contains only those symbols which are used 

throTi^hout the text. 

T Return travel time of the laser pulse. 

CO Angular velociiy of the earth rotation. 

© Greenwich apparent sidereal time, 

X, y Polar motion coordinates, y is positive westward. 

C ' Celestial pole . 

(U) = (U, V, W) Geocentric coordinate system whose body-fixed motion is 

due only to the co mm on crustal motions of the participating 
stations , 

p, ^ , A Spherical coordinates in the system (U); p Is the geocentric 

station distance; ^ and A are the station latitude and longi- 
tude. 

(U') = ([/, V', W') Geocentric coordinate system whose third axis coincides 
with the celestial pole (C) ; u' axis is along the Greenwich 
mean astronomic meridian. 

p,^'. A' Spherical coordinates in the system (U^) ; A^ = 0 is in the 

Greenwich mean astronomic meridian. 

(tf ) = {Xf't V\ w " ) Geocentric coordinate system whose third axis coincides 
with the celestial pole (C) . The U^-axis is nearly body- 
fixed . Its body- fixed position is a function of an error in 
Greenwich apparent sidereal time © and right ascension of 
the lunar refleetor as used in the computations. 
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(X) = (X,Y, Z) Geocentric coordinate system whose third axis coincides 

with the celestial pole (C) . The X-axis is along the direc- 
tion of the true vernal equinox. 

A Geocentric reflector distance. 

a, 6 Right ascension and declination of the reflector. 

[X] Set of estimable parameter combinations . 

A Design matrix for estimable parameter. 

A Design matrix in analysis of variance models. 

[Y] Set of all parameters (estimable and non-estimable) . 

A Design matrix corresponding to parameter set [Y] . 

H Matrix which transforms non-estimable parameters to 

estimable ones. 


To 


Epoch of the initial interval. 
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1. rNTRODUCTION 


Lunar laser ranging (LLR) is on the verge of -becoming- a widely 
used tool in high precision geodesy. Observations have been success- 
fully carried out at the McDonald Observatory at Austin, Texas for 
about six years. Ranging occurs regularly to the reflectors placed on 
the moon by Apollo 11, 14 and 15 and to the reflector on Lunakhod n. 

The Apollo 15 reflector is the most favorable to acquire since it is 
the largest among the lunar reflectors and it is located dn close vicin- 
ity to distinctive surface features which make easy guiding of the laser 
possible. The accuracy of ranging has been increased alreatty during 
its time of operation. A typical accuracy at this time is about 15 cm 
and better [Mulholland, 1975]. This measure is based on several re- 
turns which form a normal point [Abbot et al. , 1973] . It is basically 
the laser pulse width, the electronic calibration, and the finite size 
of the reflectors which limit the ranging accuracy. A large number of 
parameters have reportedly been improved using lunar laser ranging 
data. Because the presently available observations are being made at 
the same observatory, one cannot solve yet for all parameters which 
are of scientific interest. In the near future, stations in Hawaii, 
Australia, Japan, France and possibly West Germany are expected 
to start a lunar laser ranging program. 

Besides the stations mentioned above, which are all of the obser- 
vatory type, i. e., fixed to their respective locations, mobile laser 
stations may become" available shortly. In order to ensure the optimal 
use of range observations, it is necessary to design an optimal sta- 
tion distribution for the mobile lasers to fulfill the objective of the 
observation campaign. The analysis presented' by Silverberg’et'; al'. 
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[1976] already indicates a possible function of the transportable laser 
station as a means for densification of accurately determined geo- 
centric station pdsitions7 The ""positions of~the fixed' observatories are" 
regarded as fundamental points and the mobile stations can be set up 
in geophysically interesting areas. This procedure requires good earth 
orientation parameters (polar motion and earth rotation variation),.. 

This information is supposed to be provided by the fixed observatories. 
However, the fixed observatories are also likely to be engaged in pro- . 
grams designed to improve the models about the lunar orbital motion 
and lunar rotation (libration). Such an observational program includes 
systematic ranging to all. reflectors and is not necessarily the best- 
.program to give the earth orientation parameters frequently. It is, 
.therefore, desirable to establish an earth orientation service indepen- 
dently- of the fixed observatories. It is the subject of this stu(fy tto 
investigate the feasibilily of mobile laser stations to provide such .an 
earth orientation- service. ^ - 

. The method to be investigated is that of range-differencing where- 
by a new observable is formed by differencing the range measurements 
of two widely separated stations. These two co-observing stations • 
constitute one observational unit. As for terminology, the two sta- 
tions form the end points of a "line./' In this study, the influence 
of the length, location and orientation of a line on the recoverable 
accuracy of the orientation parameters will be investigated. . For the 
ideal case of simultaneous observations (same reflection time of the 
two pulses at the reflector) the range difference is very insensitive 
to changes in the- geocentric reflector distance because of the small 
angle subtended by the paths of the two pulses. A change in declina- 
tion, or right ascension, on the contrary, has a strong effect on the 
range difference depending on the orientation of the line. It is • 
assumed throughout this study that a good lunar ephemeris is available. 
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It may be a numerically integrated ephemeris or an analytical ephem- 
eris. Since the primary objective of this study is the design of a 
method .to determine the earth orientation and not the motion charac- 
teristic of the moon, no attempt will be made .to separate the lunar 
orbital motion from the rotational motion (libration). It is therefore 
sufficient to .carry out the range observations to only one reflector, 
say, the Apollo 15 reflector. The spatial position of the reflector will 
be parametrised by its geocentric distance, its declination and right 
ascension. The approximate values of these positional elements are 
computed from the available ephemeris and libration model. Constant 
corrections to the declination and right ascension of the reflector will 
be solved for daily. For nearly simultaneous observations, range-differ- 
encing occurs between those two observations whose reflection times 
are closest, hi such cases, the requirements on the lunar ephemeris 
are more stringent as is explained in Section 2.3.1. 

The coordinates of the stations will be transformed to their dif- 
ferences and sums. Range differences allow for the accurate determi- 
nation of coordinate differences. However, coordinate differences 
completely determine the orientation of the earth . It is pointed out 
that lunar laser ranging gives the orientation of the earth only, rela- 
tive to the lunar motion. If the results of limar laser ranging are. to 
be related to a frame defined by quasars, complementary observations 
such as differential VLBI (very long base line interferometry) obser- 
vations of the Apollo lunar surface experiment packages (ALSEPs) 
and quasars are needed. 

The earth orientation parameters will be solved for each interval. 
Since the progressive Chandler motion is smaller or equal to 10 cm 
per day, it is natural to limit the length of the interval to one day 
for which the orientation parameters can be taken as a constant. 
Actually, the length of one interval depends on the station distribution. 
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Each interval includes only one lunir transit. The interval begins 
'With-the- -epoch- of- the- -first -observation at the most‘-eastern~stat-ion- 
pair and ends with the last' observation at the most western station 
pair. The shorter the interval, the simpler is the mathematical 
model of the adjustment and the number of solution parameters is 
minimized. However, a short interval requires a favorable station 
distribution and an adequate observation schedule in order to' provide 
sufficient geometrical strength and enough observations for the param- 
eters to be deterinined accurately in each interval. 

- . • ' 

In the subsequent section, the observation equations for range- 

differencing are set up assuming that the range observations are al- 
ready corrected for influences of atmosphere and solid body earth 
tides. Fortunately, the atmospheric correction can be computed ac- 
curately to i 1 cm for the wave ler^h of the lunar lasers [Mulholland, 
1975]. The solid body earth tides can reach an amplitude for the yer- 
tical displacement of 30 to 40 cm. It will be necessary to compute 
these displacements as accurately as possible using elaborate geo- 
physical models. The complete description of the observation equa- 
tions also requires relativistic considerations. Without taking recourse 
in great detail to relativistic theory, the pertinent procedure is given 
for computing the exact time delay. The mathematical formulation 
used in' this study is such that not only chords and angles between the 
stations are estimated but also the orientation’ angles relative to an 
epoch To which is the zero point of counting. This procedure is dis- 
tinctly different from what is usually referred to as hmer Constraint 
in which case the resulting "orientation parameters" do not allow a 
simple interpretation. Due to the lack of real observations, an analy- 
sis of variance is made to investigate the capabilities of the proposed 
method as a function of the number of stations, their distribution and 
the observation schedule. The analysis of variance model is derived 
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from the rigorous model by making certain simplifications. One of 
the simplifications is to neglect the earth rotation during the travel 
time of the pulse, or equivalently, an infinitely large velocity of 
light is assumed so that the instants of transmission, reflection and 
reception are the same. In other analysis models, the terms of the' 
characteristic order 1/60, i.e. , those terms which are a function of 
the ratio of the geocentric station distance to the geocentric reflector 
distance are ne^ected, and the declination is taken constant during 
one interval. Special attention is given to the minimum number of 
stations required in order to uniquely determine the earth orientation 
and to the difficulties introduced by crustal motions. Any data loss 
due to weather conditions is not taken into account. If for a particu- 
lar interval no determination of the orientation parameters can be 
made because of loss of data due to adverse weather conditions, the 
parameters can be interpolated from their values in the adjacent in- 
tervals. 
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2. MODELLING RANGE OBSERVATIONS 


The direction of the celestial pole (C)is best suited as a refer- 
ence direction for the third axis of the coordinate system in which 
range equations are formulated, hi fact, the body motions of the ce- 
lestial pole (C) are directly estimable from laser ranges. In Leiek 
[1978], an extensive discussion on the implications of the differences 
between the instantaneous rotation axis and the celestial pole is given. 
The principal properly of the celestial pole is that its direction has 
neither body- fixed nor space-fixed periodic diurnal motions, in this 
section, the rigorous least squares formulation for range-differencing 
is given. The estimable parameters are identified,, and the earth 

I 1 •, 

orientation parameters are specified on the basis of an orthogonal 
transformation (over determined case). 

2.1 Time Delay Computation 

Processing range data requires the precise computation of time 
delay, i.e., computing the elapsed round trip time for the pulse from 
the approximate parameters. Such a computation makes relativistic 
corrections necessary. The corrections to be considered here are 
undisputed in relativity theory. Although there are several competing 
relativities theories, each gives the same corrections, at least to 
the accuracy which is needed for laser ranging within the solar 
system. The following time scales need to be distinguished: 

a) coordinate time t 

b) proper time r 

Coordinate time is the time argument in the ephemerides of the 
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planets, the earth-moon bary- center and the moon in the heliocen- 
tric coordinate system. These ephemerides are, of course, based 
on the re’ativistic equations of motion which are obtained from the 
Newtonian equations by amending relativistic perturbative accelera- 
tions. For the purpose of this discussion, one can consider the time 
scale of the atomic clock on the earth as a realization of the proper 
■ time scale. 

In relativity- the ratio of the intervals of the coordinate and . 
proper time scale is not constant. The proper time scale is a fun- 
ction of the total potential U at the position of the clock and of the 
solar system barycentric velocity s of the clock. The relativistic 
corrections for laser ranging in the solar system are developed in 
great detail in Moyer [1976] from where the following differential 
equation is taken which relates proper time and coordinate time: 

.ll = 1 _ - i/iV 

dt c® "^c/ (2.1-1) 

This is ,an approximate expression in which terms of order 1/c^ are 
ignored, c is the light velocity. Since U is positive, it is seen that 
for a fixed atomic clock on earth dr < dt, proper time r falls behind 
coordinate time. Since 

dN 

■ V 

where dN is the nurhber of observed cycles and n is. the number of cy- 
cles per second atomic time (conversion factor), one can make- the in- 
terval of proper time agree on the- average with coordinate time by 
selecting the appropriate n. This procedure was, in fact, followed 
when the TAI (International Atomic Time) second was defined. The 
remaining periodic variation between proper and coordinate time is 
derived by integrating equation (2.1-1). Moyer [1976] gives the 
following expression: 
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t ~ T = Ar 


+ 1.658 X 10 ® Sin E 

"+' 3‘. 18 ■xl0"^° pcos '4»' sm“(UTI+ ’A:) 

+ , , , ■ ■ - ( 2 . 1 - 2 )' 
The first term is the constant of integration. ' The lAU (1977) adopted 
At= 32.184 s 


It basically represents the constant shift between atomic time “TAl and 
ephemeris time ET at the initial epoch of atomic time. The second 
term is the largest periodic term. Its magnitude is about 2. ms. E 
is the eccentric anomaly of the heliocentric orbit of the earth-moon 
bary center.' The third term has a magnitude of about 2 /as. p cos ^ 
denotes the spin axis distance of the clock in kilometers, UTl stands 
for universal time, and A for the station longitude. The complete ex- 
pression contains many more terms including terms which are a ' 
function of the position of Jupiter and Saturn. 

Given the (t - T) correction of equation (2.1-2) coordinate time 
can be transformed to proper time at any epoch and vice versa. The 
computation of the travel time of the pulse is complicated by another 
relativistic phenomenon, which is called "radar time delay." This de- 
lay is a function of the light propagation characteristics as the light 
travels through the potential field of the sun. .The theoretical back- 
ground of this phenomenon is given in [Misner et al., 1973, p. 1103]. 
Computationally, the radar time delay is expressed by the so-called 
"Light -time equation (LTEq.)." Moyer [1976] gives 



rki 2GS ri, + ri + ri-^ 

+ — 3" In 

c c rj, + r^ - ri^ 


(2.1-3) 


with 

% - II?; (T,) (T,)l| 

= 11% 

ri = Ilfs (Ti) II 
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The superscript s denotes heliocentric position. GS is the heliocen- 
tric gravitational constant. The light time equation thus relates the 
coordinates of two points, k and 1, to the coordinate time t for light 
to travel from one of the points to the other. This equation has to .be 
solved by iteration. In case of LLR, the up and down legs have to 
be solved separately. 

In summary, the computed range observable Tq is arrived at 
according to the following scheme: 


Ts' 


r 

(t ' tL 




LTEq 


->-T„ 


LTEq 

>-Ti 


and r, = Tg - T^ - (t - + (t. - --tL (2.1-4) 

The subscripts 1, 2, and 3 denote the instants of transmission, re- 
flection and reception of the laser .pulse. Thus, Xg is the observed 
reception time of the pulse. It is cmverted via equation (2.1-2) to 
coordinate time Tg, which is used together with an estimate of T^ in 
order to solve the light time equation (2.1-3). The solution gives an 
improved value for Tg so that the solution of (2.1-3) can be repeated 
until no significant change in Tg occurs. The same iteration is done 
for the downward leg. The light time equation has to be solved twice 
because the upward and downward path of the laser pulse is at a dif- 
ferent position within the gravitational field of the sun. The third and 
fourth terms in equation (2.1-4) convert the round-trip light time from 
an interval of coordinate time to an interval of proper time at the 
epoch of observation. Finally, the computed time delay, , and 
observed time delay on the atomic clock, Tt - r^, can be . 
compared. 
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2.2 Basic Observation Equation 


2.2.1 Single Station Observations 

Two solutions of the light time equation (2. 1-3) and the relativistic cor- 
rections (2. 1-2) have resulted in heliocentric positions of the reflector, the 
geo center, the laser station and in the computed time delay To. This computa- 
tion is, of course, .based upon the approximate values of reflector position, 
earth orientation parameters and laser station position. On the basis of the 
observed time delay, some of the parameters, particularly those for earth 
orientation, can be improved. 

In Figure 2.1 the symbols C and R denote the celestial pole and the 
lunar reflector, respectively. The numbers 1, 2, and 3 denote the instants of 
transmission, reflection, and reception of the laser pulse. The observation 
equation will be expressed in a geocentric frame whose third axis coincides 
with the celestial pole (C). It is readily derivable from the light time equation 
as ajpplied to the up and down leg. Denoting the second term on the right-hand 
side of equation (2 .1-3), which is die contribution to the light time from gen- 
eral relativity, by Ajdthe observation equation becomes 
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0 


( 2 . 2 - 1 ) 


|Xr2-Xx|| + IIXrs-XsI! - cT - cw = 
where • 

W — - Ai 2 - -^23 + (t - ’^)tg “ (t - T)tj^ 

is a small but computable term. T is the time delay. The coordinates refer 
to the system (X), the first axis being located along the true vernal equinox. 
The light velocity c determines the scale of the configuration. By dividing 
equation (2. 2-1) by c, each coordinate is expressed in units of light 
velocity, so that any changes in the adopted light velocity result in a 
computable change in the coordinate length (scaling). 

The station position at the instants of transmission and' reception 
are related by a rotation around the third axis as follows;' 

% = Hq (COT)Xi 

where io is the angular velocity of the earth rotation. The rotational 
position of the earth is introduced by 

Xi = Ea (-©i) U' 

where © is the apparent Greenwich sidereal time at the instant of 
transmission and (U ') is the coordinate system whose third axis .still 
coincides with the celestial pole but whose first axis' is fixed to the 
Greenwich mean astronomic meridian. Substituting these equations 
in equation (2.2-1) gives the basic mathematical model: 

F(L^,XJ= ilx ,2 -Rs(-0i) !U 1 !Xr 3 -R3(C0T -©i) ifll - cr -cw = 0 

( 2 . 2 - 2 ) 

This is the standard model in adjustment theory for the case where 
observed quantities and unknown parameters are in the same equation. 
With the notation of Uotila [1967], the linearized form of. (2.2-2) is 
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(2.2-3) 


B V + A X + ,W, = 0 

r n r u . i 

where 

B = —■ , A = “, andW = F(L^, Xo) 

O La ° 

r denotes the number of equations and u is the number of parameters. 
The weight matrix is the inverse of the variance-covariance matrix 
of observations multiplied by the variance of unit weight, 

P = CT ® Zj 
0 Lb 

As usual, the subscripts a, b, and p refer to adjusted, observed and 
approximate values, respectively. The partial derivatives in B and A 
are evaluated for the approximate and observed values. V denotes the 
residuals. X is the set of solution parameters; they are corrections 
to the approximate values . hi order to avoid confusion, the parameters 
will sometimes be denoted by [X] . The usual minimization of 
y' PV 

gives the least squares estimate for the parameters [X] 

X = - (A' M'^ A) A ^ M’^W ' (2 . 2-4). 

with 

M = B P'^ 

The adjusted variance of unit weight is 

-2 _ V'PV 

- " r - n (2.2-5) 

Finally, the variance-covariance matrix of the adjusted parameter is 

S. = 6-^ (A^M'^-A)’^ (2.2-6) 

The partial derivatives of the function F will be given below for 
spherical and Cartesian coordinates. Both systems are, of course, 
strictly related, but in some cases of analysis of variance, one or 
the other system is more convenient to interpret. 

The subscript "2" on the lunar positional elements is subsequently 
omitted. Unless otherwise stated, its position is always evaluated at 
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reflection time. Similarly, the subscript "1” on the symbol B is de- 
leted, remembering that the apparent Greenwich sidereal time has to 
be evaluated at transmission time. Denoting the spherical coordinates 
of the station in the (U^) system by (p. A' , the partial derivatives 
in the design matrix A are computed from equation (2,2-2) as follows: 

A ~ {A- pcos cos 5 cos (A' + 0 - a) - psin sin 6 } 

La P‘12 


•23 


{A - fjcos cos 6 cos(A^- orr + 0 - o;) - psin sin 6 } 


A = 
a 


pA 

3 

P A 


p4a ^ » X / ^ 

COS $ COS 5 sin (A +0 - a) 

ri3 ' 


■3S 


cos cos 6 sin ( - WT + © - oJ 


^6 = 


P A 
r 


— {cos sin fi cos(a' + 0 - 6 ^ -sin cos si 
X 2 - ■ ''f 


Y"— 6 cos(/^' - £jor + ©“O') “ sin cos 5 


‘'S3 


}>(2.2-7) 


A 


P 



{ p “ Acos cos fi cos(A^ + 0 - p)| 

+ { P - Acos cos fi cos (A' - 6t)T + 0 - O')! 

I sin cos 6 cos (A' + © - o) “ sin 5 cos [ 

^12 ^ 

P A f • / i 

+ __ ( sm $ cos 5 cos (A - toT +0 - o) “ sin 6 cos # | 


A / — 
A +0 



cos cos 6 sin (A^ + © - 0 ) 


+ j — cos cos 6 sin (A^ - cor+ © ~ 0 !) 

32 ^ 

Note that A' + © ~ & is the hour angle of reflector . The' symbol r 
stands for the distance between' the observation station and reflector, 
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in particular 

rs3= ||j^R2-^3ll 

It is seen that the coefficients A and A. / _ are linearly dependent. 

Cl A 

Thus, only the linear combination + © - o:, i.e., the hour angle of 
the reflector, is estimable so that the list of parameters contains five 
estimable quantities 

[X] = [dA,d5, dp, d>^', d(A' + © - a) ]. (2.2-8) 

Note that parameters or linear combinations of them are called estim- 
able if the corresponding design matrix A is non-singular. 

The Matrix B is diagonal since each equation of (2.2-3) contains 
only one observation. The diagonal term is 

B = cos cos 6 sin (A'-cor+©-o')-c 

r r33 

Each observation from the same station adds one equation to the 
system (2.2-3) and three parameters, concerning the reflector position 
and the earth rotation, to the list of (2.2-8). For k observations, the 
complete parameter list is 

pq = [dA\ d6^ ... dL^, d5^ dp, d^« , d^(A' + ©-a), d®(A' + 0 - 
... d*'(A' + 0 - a) ] 

A least squares solution becomes possible by imposing the constraints 
d/i’^ = dA 


d6^ = d5 

dcP = do; 


(2.2-9) 


dd^= d0 

for r) = 1. ..k. These constraints make it necessary that a good lunar 
ephemeris and station clock are available in order for such simple 
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modelling to be permissible. A detailed discussion on these constraints 
follows in Section 2.3.1. If one incorporates the constraints, there 
remain five parameters which can be solved by least squares 

[X] = [d A, d 5, dp, d d(A' + 0 - a ) ] . (2.2-10) 

The station positions in (2.2-10) can be interpreted as referring 

to a system (U^^) which differs from (U^) by a small rotation around 
the third axis. The rotation is due to errors in the clock, d ©, and 
lunar right ascension da. Since .the station coordinates in the (uVsystem 
are a function of time due to polar motion, a reformulation qf the mathe- 
matical model is given in terms of station coordinates in a conven- 

tional geocentric terrestrial system (U), which coincides with the (U^)- 
system at some standard epoch, 

(U) = (U'') at epoch To . 

The relation between these systems is given by the polar motion co- 
ordinates (x, y). The same representation as found in [Mueller, 1967, 
p. 82] is chosen; i. e. , the origin of the polar motion coordinate sys- 
tem is at the pole of the (U)-system, the x-axis is, .along the direc- 
tion of the zero longitude, and y points westward. The latitude and 
longitude in the (U)-system, ^ and A, are related to those in the 
(UVsystem by [Mueller, 1967, p. 87] 

d((|»- = y sinA - X cos A 


d(A- A^) - -(X sin A + y cos A) tan ^ 

The partial derivatives for latitude and longitude in the observation 
equation (2.2-3) now become 


A /d® + Aw 


9 


A'+e-a 


dA 


= A-/ [d^ - y sin A -i- x cos A] + A. / ^ [ d(A O') + 

9 A •?'© -Oi 

+ (x sinA + y cos A) tan 
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( 2 . 2 - 11 ) 


= Aw - Aw y sin A+ A . / x cos A+ A , A 
<p |) A +0-0! 

+ A ./^ ^ X sinA tan^+ A _ y cos Atan ^ 

A -H9- a A+B-cx 

The coefficients (2.2-7) are evaluated with and A' , i.e., with the 
approximate station coordinates in the (U')-system. If those latitudes 
and longitudes are replaced by $ and A then one formally obtains the 
partials which one would have obtained if the station position in equa- 
tions (2,2-2) had been expressed in the (U)-system, and the partials 
had been evaluated with x = y = 0. No further changes are needed in 
the coefficients (2.2-7). 

The parameters of equation (2,2-11) together with those of (2.2-8) 
give the following list of estimable parameter combinations: 

[X] = [dA, d5, dp, p, V] (2.2-12) 

with 

fX = d. = d^-y sin \ + x cos A 

l^=d(A^+ pi ~ <y) = d(A+©“ a) + x sin A tan $ + y cos A tan 


The coefficients of the last two parameters are 

A = A^ 

A, I A . — 

^ A+0-a 

There exist an alternative set of estimable parameters. From 
equation (2.2-11) the coefficients for polar motion are 


A = A , cos A + A , ^ sin A tan $ 

X 9 A+©-a 


A = -A,sinA+ A, * cos A tan ^ 
y 9 A+Q-a 


(2.2-13) 


Incorporating the inverse relations of (2.2-13), 
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= Ax cos A - K sin A 

$ 

^A+0- a "" ^ ^ ^ ^ 


in the general observation equation gives the followir^ set of estimable 
parameters: 

- [dA, d6, dp,ij,\ v'] (2.2-14) 


with 


and 


ji' =x + cosAd*J)+ eot'S'sin A d(A+ © “ a) 
~ y - sin Ad^ + cot $ cos A d(A+ 0- ^ ) 

A / “Ax A,,/ — Ay 


It is emphasized again that all coefficients have to be evaluated with 
the approximate station coordinates in the (U) -system. 

Analogous expressions to (2,2-7) for Cartesian station coordinates 

are: 


A. - { -U^ cos 5 cos (©-a) + V'' cos 6 sin (©- a) -W^sin6+ A] 

^ r^s 

C~U cos 6 cos(-uj r+ © - 0 ! ) + cos6 sin(-6(jr+ 0- Q:) 

^33 

- W'^ sin 6 + A } 


A - Cu^ Sind cos (0- a) - sin 6 sin (0- a) - cos 6 } 
0 ^10 

+ — {U^ sin 6 cos (-m t + 0 - a) - sin 6 sin(-cor+ 0- a) 

Tag 

- w' cos 6} 
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^ cos 6 sin (e~ od “ ^ 

ri2 


A, ^ — 
a-0 

+ A (-TT^ COS 6 sin(-ooT + 0 - o:) -v" cos 6 cos(-wt+ ©-«)} 

. I'ss 

A , = — fu' - AcosS cos (0 - a)} -i- J_{u^ - Acos6 cos(-^*)T +Q-a)} 

ri3 ^ raa- 

Ay/ = — tv' + A cos 5 sin (©- a)} i— [v' + Acos 6 sin (-wr + 0 -Q:)} 

i"ig r33 

t 

A.-/ = — — tw'- Asin 6 } + - AsinS } 

3^13 ^33 

The diagonal element of the B- matrix is 

B = Ah f-u' cos6 sin(- 60 T + © -O) - y' cos 6 cos(-iOT + © - a)} - 
T ' • 

Combining the linearly dependent coefficients A^_ 0 jAii/ , and A-.-/ 

A = -u' Ay/ + V’ Au' 

O'- 0 

one finds again a list of five estimable quantities 

[X] = [XiIXs] = [dA, d 6 ; du'\ dV' , dw" ] (2.2-16) 

with 

d u" = d U' + V ' d(Q! - 0 ) 

d V" = d V' - ,U' d(a: - 0) 

dw" = dw' 
and the coefficients 

Ai;// “ Ay^ } — A'y/ 

The constraints (2.2-9) have been incorporated in the list (2.2-16). 
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The expressions are transformed to those of the (U) -system by 
introducing polar motion coordinates as follows (only terms of first 
order are retained): 

U' U - xW 

V'=V + yW (2.2-17) 

W' = xU - yV + W 

Differentiating equations ' (2 . 2-17) and combining them with the param- 
eters (2.2-16) gives another set of estimable quantites: 

[X] = [XiiXs] - fdA, d5i dt/', dv'', dw"] (2.2-18) 



The coefficients are 

Aq''' = Ay » A',// = A, » Af^/ = Ai^ 

Note that the coordinates (U^ of the coefficients (2.2-15) have to be 
replaced by (tO • The partials for polar motion, Aj^, and Ay, are not 
needed explicitly although they can be derived" in a similar manner as 
was done for the case of spherical coordinates (2.2-11). For reasons 
of abbreviation, the rotation parameter in (2.2-19) is denoted by a, thus 

do: = d(o: - 0) 

The quantities dU'^, dV'^ and d can be interpreted as coordinate cor- 
rections in the estimable frame of reference. The third axis of this 
frame coincides with the direction of the celestial pole (C) , whereas the 
first axis deviates from that of the (U ^system by an angle da. 
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Equation {2.2-19) can equivalently be written as 

•d#' ■= -dU + [Es (d a) (y) Es (x) - I ] U 

The derivatives so far have resulted in a set of estimable param 
eters for laser observations at one station. The columns of the de- 
sign matrix A were checked for linear dependencies and the corres- 
ponding parameters combined in order to yield estimable parameters. 
This concept is identical to that of rank factorization. For the sake 
of subsequent discussions, the identities between both approaches are 
pointed out. Let [Y] be the vector of all the m parameters which 
have entered the formulation 

[Y] = [XilYglYg] = [dA, dd: x , y, da-dU, dV, dW] 

The corresponding singular design matrix A of size (r x m) is 

rA„ = (Ai i Aa i As) 

The rank is 

R(A ) = m - s = u 


with rank deficiency s = 3 for the present case. [Yg], which denotes 
the parameters of the earth orientation, still contains the combina- 
tion of right ascension and Greenwich sidereal time. Both could, of 
course, also be separated, and thus increase the rank deficiency to 
four. In linear algebra it is proven that the matrix A can be fac- 
torized such that 


and 


rA„ = 5 .D^^ uHn with r >u, m > u 


Thus, 


R(D) = u 


AY = DX 
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gives the estimable parameters [X]. The following three identifica- 
tions can be made: 

1) The set of estimable parameters is the same as in (2.2-18), 

2) The design matrix for the estimable parameters is D = (Ai-rAs), 

3) The non-estimable parameters are transformed to estimable 


quantities by the 

matrix 





r- 

1 

0 

0 

0 

0 

0 

0 

0“ 


0 

1 

0 

0 

0 

0 

0 

0 

H = 

0 

0 

-w 

0 

V 

1 

0 

0- 


0 

0 

0 

W 

-u 

0 

1 

0 


.0 

0 

u 

-V 

0 

0 

0 

1_ 


The H matrix is written for abbreviation as 


Evaluating 


H = 


TOO 
O 1 


A = DH 


— (Ai A3) 


I O 
O F^ 


0 

1 


— (Aj^ : A3 F A3 ) 


yields the linear relation between the coefficients; 



2 . 2.2 


Range Differencing 


The adjustment model for range-differencing is readily derived 
from the basic model (2.2-3) and its set of estimable parameters of 
equation (2.2-18). Denoting the two co-observhqg stations by i and j, 
the linearized form of the adjustment model is 


BtVi + A^Xj + Wi = 0 


( 2 . 2 - 20 ) 


BjVj + AjXj + Wj = 0 

Since these equations included the constraints (2.2-9), there are 2r 
equations for u = 8 estimable parameters, i.e., the parameters 
dA, d6, dck: are common to both sets. Subtracting the first equation 
from the second in (2.2-20) gives the new set of equations. 



■vr 


-X," 

(-Bi Bj) 


+ (-AiAj) 



-Vj- 




( 2 . 2 - 21 ) 


with 

Wij = Wj - Wi . . 

It is understood that in (2.2-21) those terms are differenced 
which correspond to the observation equations for simultaneous or near 
simultaneous range observations at the stations i and j. This is 
automatically achieved if the equations in (2.2-20) are ordered suc- 
cessively in time. The parameters are conveniently transformed 
to their sums and differences as follows: 


(-Ai Aj) 


X, 


X, 


(AMj) • 


-Xi 




+ (Aj - Aj) — 1 


X, ' 

2 


With the notation 


- X, -X, ^ _ X, + X, 

’ Jt* 2 


( 2 . 2 - 22 ) 
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the equation (2.2-21) now becomes 


(-BiBj) 


Vi 


+ (Aj + A, ! A, - Ai) 


‘3-i 


X 




+ Wij = 0 


The list of estimable parameters in the above equation is 
[X] = [Xj Xo] . 


= [dA, d6 i djfyt, dVj^,, dl/;+., dVUh dW^+i] ' 

(2.2-23) 

and the design matrix consists of the differences and sums of the 
original coefficients, 

(A) = (All A 3 ) 


(Aj^ Aj^, Aj 6 -A.i5 


Aju + A 




Ajv A(v» A 1111 A 


^iVj 


Jw 


^iw> 


Aju - Aiu> Ajv - Aiv, Aji^ - Ajw) 

(2.2-24) 

Note that the geocentric reflector distance and the declination in 
(2.2-24) are not split up into their sums and differences because 

t ■ 

dA j - dAi = 0 and ‘d 6j - d 6 1 = 0 is valid throughout the interval 
according to constraints (2.2-9). The last six parameters in (2.2-23) 
are the corrections to the sums and differences of the station 
coordinates in. the (U")-sy^em. They are. related to the corres- 
ponding parameters in the conventional terrestrial system (U) by 
the three orientation angles x, y and da . With equations 
(2.2-19) and (2.2-22) one gets the following expressions: 


1 

p 


du^i 


■* 

0 

da 

-X 


Uj-i' 

-dv';_i 

1 

■ =; 

dVj-t 

+ 

-da 

0 

y 


Vj_i 



dWj,-.. 


X 

-y 

0 


Wj-t 


(2.2-25)--' ■ 
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and 



{2v2-26)- - 


To the first order of the orientation parameters these equations can 
be written as: 

dU = dUj_i + [Eg (do; ) Ri {y)Rs (x) - 

du = dUjfi + [Rs (da) Ri (y) Rs (x) - I] 

I is the identity matrix. The least squares estimates according to the 
general equations (2.2-4) to (2.2-6) are, with 


M - Bj + Bjl^bj Bj), 


4 ^ 

T T 

Ai^M’^Ai 

AlM'^-Ag 

-1 

A^M'^ 

[Xf= [Xl, Xgj" = 

_A3M% 

Al M-^Ag_ 


Ag" 


Taking the inverse analytically yields 
■ ‘ Xi = GiWij (2.2-27) 

where 

Gi =■ --{a 1 Ai- Ai^ M-'^Ag (Ag^ M%) A^g M'^Ai 3 Aj M'^ 

+ (A\ M'^Ai ) ‘^A J M '^Ag Q X Ag M 

3 


Xa = .Qx^[A^M-"Ai(AIm-'Ai)-'aIM'^ - A^ M‘'3Wt3 (2.2-28) 

with 

Qx- = IM M'^Ag - Al M''Ai(Al M'^Ai)'" Al M'^Aa}'^ 

3 

and • , 

3 

Sx = %Qxg (2.2-29) 

3 

whereby ao is computed from (2.2-5). 
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As for rank factorization, the necessary identifications are 
readily available. Let [Y] include all the parameters 

[Y] = [Xi i Y2 i Ya ] 

= [dA, d6 : x’, y,da : dllj-^, dWj_t, 

( 2 . 2 - 30 ) 

dV^i, dW^i ] 

with design matrix 

A - (Ai : Ag : A3) 

which has a rank deficiency of three, then the rank factorization is 


AY = DHY = DX 
X = HY. 


The design matrix D for the estimable parameters is 

D = (Aj : A3), 

where 


R(D) = u. 


The non-estimable parameters are transformed to a set of estimable 
parameters by 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

?. 

0 

0 

0 

0 


0 

0 

-Wj_i 

0 

Vj_i 

1 

0 

0 

0 

0 

0 

0 

0 

0 

Wj-i 


0 

1 

0 

0 

0 

0 

0 

0 

U 3-1 

-Vj-i 

0 

0 

0 

1 

0 

0 

0 

0 

0 


0 


0 

0 

0 

1 

0 

0 

0 

0 

0 

Wj-H 

-UiM 

0 

0 

0 

0 

1 

0 


0 


-Vifi 

0 

0 

0 

0 

0 

0 

1 


For abbreviation H is written as 



O O 
F^ I 
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( 2 . 2 - 31 ) 



The estimable parameters [X] are those of equation (2.2-23). Further, 


A- = DH 


- (Ai As) 


I 

O 


o o~ 


— (Ai : AsF^ : A3) 


yield again a relation between the coefficients 

As = A3F^ 


2.3 Aspects of Implementation 

2.3.1 Ephemeris Modelling 

The adjustment formulation includes the constraints (2.2-9). 

The first three constraints imply that for the length of an interval 
the corrections to the reflector position are modelled by three con- 
stant parameters. Such a constraint makes simultaneous obserya- 
tions superfluous, at first sight. It requires that a good lunar ephem- 
eris be used to compute the reflector positions at the instant of refleo* 
tion of the laser pulse. This simplified modelling seems justified in 
view of the smooth motion of the 'moon due to its large moment of 
inertia. The shortest libration terrn is 13.6 days. Even the presence 
of short periodic terms does not invalidate the above constraints as 
long as the coefficients of these terms are correct. Since the lunar 
orbit is quite accurately derived from rigid body theory, the author is 
not aware of any frequencies, say, analogous to the critical frequen- 
cies in the nutations of the earth as a result of core motion, that are 
significantly wrong in the available conventional ephemerides. 

The constraints appear even more reasonable if one accounts 
for the fact that a given range accuracy is capable of resolving the 
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reflector position to only a certain level. Figure 2.2 displays the 
geometry of the reflector and stations i and j. The earth rotation 
during the travel time of the pulse is neglected for simplicity. 
Figure 2.2 is valid for either a north-south or an east-west line. 



Figure 2.2 Resolution of Reflector Position 
The expressions for the topoeentric distances 

r® = p® + - 2pAcos 4 

r® = fP + ^ - 2p A cos (17 - 4) 

can be expanded in terms of in order to give the range difference 
as follows: 

Tj --ri_= p [cos i - cos (77 - ^ ] 

Differentiating this equation with respect to the lunar position ^ 
results in the error estimate 

[d(rj - n)| ^ 2 pd ^ 

In terms' of the linear distance at the moon, this estimate becomes 

j d(rj - rj)| ^ ^ 

oU _ Ltn 

If a range difference accuracy of 3^/~2 cm is assumed, then the 
linear reflector position can be determined at best with an accuracy 
of ±1.3 m. In actual computation, several parameters are solved 
simultaneously. The existing correlations between the parameters 
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tend to increase the uncertainty even more as compared to the simple 
calculation above. From this point of view, the modelling of dfi and 
’da as a constant per interval seems adequate. The modelling of the ' 
geocentric reflector range error d A is more critical if strictly simul- 
taneous observations are not possible. For near simultaneous obser- 
vations, the condition dAj = dAi = d A requires that the change in 
geocentric reflector distance between the two reflection epochs is 
computable at least with measurement accuracy since the computed 
range difference observable will directly depend on such an error. For- 
tunately this error can be decreased, theoretically, even completely 
eliminated, by scheduling the observations "as simultaneous as prac- 
tically possible. " If the analysis of actual observations indicates that 
the modelling of the reflector position in the manner described here 
is not sufficient, one can attempt to create artificial simultaneous 
observations by interpolating sii^le station ranges at certain epochs. 

In fact, the simplified interpolation method which leads to the construc- 
tion of the "normal points" [Abbot et al. 1973] might still be sufficiently 
accurate. In any ease, the investigation on the proper interpolation 
method should be carried out with real data and not with simulated 
observations. 

2.3.2 Timing 

The condition d^= d © an (2.2-9) expresses a perfect synchroni- 
zation between the two co-observing station clocks. Through frequent 
comparison with transportable clocks, this condition can be fulfilled 
quite accurately. If both stations are capable of utilizing LORAN C 
transmissions, it is possible to maintain a long term clock synchroni- 
zation of lys (ground waves). This corresponds to an equatorial 
rotational motion of the earth of less than a millimeter. Time syn- 
chronization errors are, therefore, virtually negligible and the 
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parameter ia, indeed, contains only irregularities of the earth rota- 
tion and errors of the Ittnar right ascension. 

2.3.3 Coordinate System Definition 

The estimable quantities of (2.2-23) are the corrections to co- 
ordinates in the (U^') -system, whose third axis coincides with the 
celestial pole (C), and whose X-axis differs from the Greenwich mean 
astronomical meridian by da. Range observations give coordinates, 
whereas for range differences, the parameter set is preferably trans- 
formed to coordinate differences and sums. The origin of the (U^V 
system is at the instantaneous center of mass. This is so because 
the earth rotates around its center of mass. Both requirements, 
i.e., origin at the instantaneous center of mass and alignment of the 
third axis with the celestial pole, are operationally achieved by ex- 
pressing the station positions as follows: 


X 


cos cos 

Y 

= Ra (-©)P 

cos sin a/ 

Z 


sin 


where p is the geocentric distance and are the latitude and 

longitude of the station in the (U ^system. This formulation was used 
when setting up the mathematical model. 

The problem of measuring the orientation in space can be regarded 
as solved as soon as the estimable parameters of the participating 
observatories become available, say, in the form of a table to be 
issued every day. The reference direction is the "mean" instantaneous 
north celestial pole for that particular interval; the word "mean" refers 
to the Chandler motion of approximately 10 cm or less per day. An 
equivalent way of representing the orientation of the earth is to issue 
daily a list of polar motion coordinates of the terrestrial position of 
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the celestial pole relative to a conventional. pole. This list should also 
contain the rotation parameter do:. It is emphasized that both methods 
of representation- are-str-ictly -equivalent-. -Both -are- merely- related -by 
an orthogonal transformation (3 rotations) which leaves the angles 
and distances between stations invariant. The adjusted differences and 
sums of the coordinates in the (U^Vsystem and their variance-covariance 
matrix serve as input for a ’’second" adjustment (transformation) which 
determines the orientation parameters 

[Yg] = [ X, y, da] (2.3-1) 

The "observations" are according to (2.2-23) 

Lb =Xs = dVjLt, dW^i, dTj'Uu dWj^U 

s 

The variance-covariance matrix of observations is given by equation 
(2.2-29) 



We note that this covariance matrix is a submatrix:, of S]; x 'j it is gen- 

X 3 

erally a full matrix. The mathematical model for the second adjust- 
ment is readily given by equations (2.2-25) and (2.2-26) which, with 
the help of submatrix F in (2 .'2-31) and equation (2.3-1) is written as: 

Ya = -F' Yg (2.3-2) 

This is the linearized form of the , adjustment model with observation 
equations. The residuals are [Y3], i.e. , the corrections to the 
coordinate differences and sums in the (U)-system. The least squares 
estimate of the parameters [Yg] are obtained by minimizing 

• Yj Ya 
3 

Note that the minimization is based on a full weight matrix. . The 
least squares estimate is given by , the standard formula. Substituting 
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expression (2,2-28) for X 3 , one obtains 
Y3 = G 2 W„ 

= [F Qx F[AJ M'"- Ai(A"i Ai)'^AJ M'^ - AJ (2.3-3) 

3 

A. 

The variance-covariance matrix of the adjusted parameters [Yg J is 

Sy =^02 (F2J"-F")'" 
a 3 

the adjusted variance of unit weight being ' 

4 = Yl ^ Y 3 . 

DF denotes the degree of freedom. It is a function of the number of 
participating stations. The residuals, [Y 3 ], which are the corrections 
to the coordinates in the (U)- system, are 


= -FT Yg + X 3 (2.3-4) 

= {-FT [F Q ~3 FTf^F + (AT m'^Ai(ATm"Ai)"'atm”'- AJmVi j 

The variance-covariance matrix of the adjusted residuals is given by 
the standard formula in least squ^es 

Sy = “ FT Syg F 

33 

The implied condition of the second least squares solution is 

FSI3Y3 = 0 

Combining (2.2-27), (2.2-28), (2.3-3), and ( 2 . 3-4), the least squares 
estimates of all parameters are 
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Xi' 


Gi“ 



-t 

0 

A 

Ys- 

. -= ■ 

69- 

w-ij =■■ 

-0 

-[-f-S'x- fVf ' -- 





3 3 

1 


G3 


0 





X 3 


(2.3-5) ■ 


K is a characteristic of the solution (2.3- 5) that [X^] does not depend 
on any m inimi zation which occurs on [Y3]. The effects of the off-diagonal 
terms in L* on the adjusted parameters, [ Y2] and [Y3], depend on the 
magnitude of the correlations between the observations. Magness and 
McGuire [1962] derived the following limits: 


^Bln ^ue ^ ^ S uc (2.3-6) 

where L uo is the variance-covariance matrix of the adjusted parameters 
[YgJ if only the diagonal elements (uncorrelated observations) 

Puc = (S- ii)'" (2.3-7) 

*3 

are used. X ajn and Xcax iii® minimum and maximum eigenvalues 
of the correlation matrix , S, of the observation noise, 

S = p|s^^3pt 


The procedure discussed above yields polar motion coordinates, 
i. e. , the daily mean position of the celestial pole (C) with respect 
to the pole of the conventional terrestrial system (U) whose position 
coincided at the zero epoch, 3^, with the celestial pole. The defini- 
tion of the system (U) depends initially on the coordinate differences 
in the (U )-csystemof the participating stations during the interval 
Subsequently only crustal motion of the defining stations can change 
the terrestrial position of the coordinate axis. Since the residuals 
of the second adjustment, [Y3], are the coordinate corrections in the 
(U)-system, it is possible to monitor station motions by analyzing 
the residuals over a longer period of time. It is very important to 
note that common station motions will be absorbed in the orientation 
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parameters (polar motion and earth rotation) . This property is implied in the 
set-up of the second adjustment. 

The station coordinates to be used for evaluating the partial derivatives 
in the first adjustment of each interval are always the same adopted- coordi- 
nates of the system (U). Equations (2. 3-3) and (2. 3-4) express. an orthogonal 
transformation between the systems (U) and (U^^), whereby the covariance 
matrix of the "observed" coordinates in the system (U) is taken as zero 
(adopted coordinates without error). This is a special case of die more gen- 
eral transformation where a non- zero covariance matrix is assigned to each 
set of coordinates. One could, of course, consider assigning the covariance 
matrix of the adjusted coordinates at interval To to ihe adopted set which 
defines the system (U). In that case the coordinates in both systems (U) and 
(U^^) would receive residuals (corrections). However, in the event of crustal 
motion the adopted covariance matrix would become increasingly inaccurate 
leading to distorted residuals. It is, therefore, reasonable to proceed with a 
zero covariance matrix for the coordinates in the (U)-system, and thus con- 
sider the adopted coordinates as uncorrelated and haying no error, and perform 

the computations as discussed above. Strictly speaking one should interpret 

^ // 

Ys in this case as a correction to the observed coordinates in the system (U ). 

A 

But such an interpretation does not diminish the value of Ya as an indicator for 
crustal motions . 
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3. ANALYSIS OF VARIANCE 


A numerical study is carried out in order to investigate the 
expected accuracy obtainable from range-diKerencing, The basic 
time span is one interval (one day). The study is a classical covari- 
ance analysis using hypothetical observations to increment the normal 
matrix. 

The normal matrix is formed for a simplified model (analysis 
of variance model). Arnold [197.4] also discusses simplification leading 
to a somewhat different analysis of variance model. Numerical inves- 
tigations were previously reported by Fajemirokun [1971], Kaula [1973], 
and Stolz and Larden [1977], All three studies are concerned with 
single ranges. Their assumptions and parametrization vary widely. 

The first two investigations carry out the adjustment over a longer 
period of time. Kaula models polar motion by four frequencies having 
periods between ha l f a month and one month. Stolz and Larden assume 
a perfectly known lunar ephemeris. They find that the orientation param- 
eters are usually obtainable to better than measurement accuracy if 
the averaging interval is two days. 

3.1 Model Simplification 

As a first approximation, the rotation of the earth which occurs 
during the travel time of the pulse is neglected. The basic model (2.2-2) 
sinipltlied as 

r= |[XR-Xi| (3.1-1) 

The simple adjustment model for observation equation 
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La = F(Xa) 


can be used. The linearized form is 

V = AX + L ‘ 

V denotes again the residuals. L is the difference between the com- 
puted and "observed" time delay. In the. present case, only the i^T. 
verse of the normal matrix 

N"’- = <A^PA)'^ (3.1-2) 


is analyzed due to lack of real observations. P is the weight matrix 
of the observations. The coefficients of the design matrix A are ob- 
tained from those of the rigorous model, (2. 2-7), by setting 

ri 2 = rgg = r 


deleting iorin the arguments, and by dividing the resultant coefficient 
by 2. With these approximations, the coefficients of the estimable 
parameters of (2.2-12) are (spherical station coordinates): 


A.= — {A - pcos$cos6 cos(A+ 0 -a) - psin$ sin6 } 




Ag = {cos^ sin 6 cos (A+ 0.- a) -sini>cos6} 
Ap= -^Cp - Acos$ cos 6 cos (A+ © - O') } 


(3,1-3) 


A^= p A { cos 6 sin^cos(A+ © - O!) -sin 6 cos } 


A+© -a 


- PA 


cos 6 cos sin(A + • © - a) 




The polar motion coefficients are according to equation (2.2-13): 
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Ax = {sin<i>cos6 cos( 0 - a) -cos cos A sin 6 } 


Ap 

L — T | 


{ sin <l> cos 6 sin( © - a ) + cos ^ sin A sin 6 } 


(3.1-4) 


For the case of Cartesian station coordinates the coefficients of 
the estimable parameter (2.2-18) are 




^ 1 

A {-U cos 6 cos( © - a) + V cos 6 sin(© - a) -W sin 6 + A } 


Ag = — { U sin 6 cos(©- Q) -V sinQ sin (© - a) -Wcos6] 


Au = ~(u - Acos6 cos (©-a) } 


) (3.1-5) 


1 I 

Av = — iV + Acos6 sin( © - oi) } 


A„ = — [W - Asinfi } 

r , . j 

In the case of range- differencing, the coefficients of the parameters (2.2-23) 
are, according to (2,2-24) and (3.1-5), as follows; 


A^- — {-Uj_t cos 6 cos(0-a) + Vj_i cos 6 sin(©-o:) - Wj_i sinfi} 
= order 1/60 


) (3.1-6) 


2A 


{Uj-i sin6 eos(©-a) -Vj_j sin6 sin(©-o:) - cosg } J 




Au = — fUju + Acos 6 cos(© - a) }. 


r* 


2 V 

S.V . = — + Acos 6 sin(©-a) } ) 


Aw = - Asin6 } 
ri ^ 


(3.1-7) 
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'a' 2 

Au = = order 1/60 

fri ^ . 

Av = = order 1/60 

jfi r J ' 


A ^ = “W. i - order 1/60 


> 

J 


(3.1-8) 


The expressions above are valid for simultaneous observations. The coeffi- 
cients of the station coordinates are of either of two characteristic magni- 
tudes, i.e., order 1 or order l/60, because of the earth-moon geometry for 
which r/A ~ 1 and p/A ~ l/60. The coefficients to the coordinate sums 
belong to the latter group as well as . The respective parameters cannot 
be determined accurately from range difference observables. The coordinate 
difference along the third axis is also determined very weakly because its 
coefficientSw is of order 1/60 for 6=0. Moreover, Aw is independent 

j” 1 J * 

of the hour angle of the moon so that its change in magnitude during one inter- 
val is very small. Depending on the location of the two stations, there will 
be a more or less strong correlation between d6 and dw"-i. Therefore, range 
difference observables are capable of determinii^ only two parameters accu- 
rately during the time of one interval; they are according to equation (2.2-25) 


dUjli = -Wj-iX + Vj-i da + dUj-i 
dVj-i = Wj-iy - Uj^i da + dVj-i 


(3 . 1-9) 


There is no disadvantage in having some coefficients only, weakly determined. 
Their significance is found by statistical testing, and, if insignificant, they 
are deleted altogether or their approximate values can be weighted. In the 
latter case, there will be no numerical problems when inverting the normal 
matrix. Good approximate station coordinates are easily available. The 
accuracy of Doppler positioning^ say ± 1 m, will be shown to be sufficient. 
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The coefficients (3, 1-8) change their order of magnitude if the obser- 
vations are not simultaneous. For example, 

= — (Uj - Aj cosSj cos(©j - ) } 

3+1 

-~ [Uj - At eos6i cos(©i - ) } 

Using the approximations 


rj ^rt = r 

e: 6i = 6 


(Xj s Oil = a 

this coefficient becomes 


- Acos6 [cos(©j - O') -cos(®t - o;) ] } 

3+1 ^ 


sf -2 cos 6 sin 


©< + © 


^ - O' 


. ©1 - © 


sin 




m ■ 


The first term exceeds the order of 1/60 if the difference in the epochs 
of observations is 

0j - ©t >4 min , 

The same limit is found for the coefficients of V. , and W-.in (3.1-8). 
Therefore, in order to keep the coefficients (3.1-8) at order 1/60 or 
less, the range observations, which form the range difference observable, 
should occur within 4 minutes. 

Another analysis of variance model is arrived at if we neglect 
terms of the order 1/60 and assume a constant declination during one 
interval . Such a model allows us to study the effect of changes in 
declination on the parameter separation, hi case of spherical station 
coordinates, the coefficients (3.1-3) are 
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1 >> 
Ag = p .cos^ sine cos(A+ © - a) - p sin^ cos 6 


/■V-/ 

Ap= -cos ^ cos 6 cos(A+ 0- a) 

~ P sin $ cos 6 cos(A+ 0 - a) - psin6 cos^ 


A 


A+© - a 


p cos 6 cos ^ sin{A+ © - a ) 


) (3.1-10) 




The coefficients for polar motion are, according to equations (3.1-4), 

Ax = p sin cos 6 cos (0 - a) - cos St*. cos A sin 6 

(3.1-11) 

Ay - p sin ^ cos 6 sin (© - Q!) + cos $ sin A sin 5 

In order to make interpretation easier, the polar motion coordinates 
(x, y) are transformed into along-meridian and across-meridian 
components (x^, y') by 


X ^ = X cos A - y sin A 
y' = xsinA+ ycosA 

Their respective coefficients are derived with the help of (3.1-11) as 
Ayf = p sin^ cos 6 cos(A +0 - a) - p cos ^ sin 6 

(3.1-12) 


Ay/ = psin^cos6 sin(A+ 0-0!) 


Since 6 is assumed constant, equations (3.1-10) and (3.1-12) make some 
additional linear combinations possible. The part of the observation 
equation which pertains to the geocentric distance p and to the latitude 
^ can be rewritten so as to contain the spin axis distance, p cos^, 
explicitly. Using the expressions of (3.1-10), one gets 

Ap dp + A^ d^ == Asa d(p cos 4>) - p sin 6 cos ^ d^ 
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where the coefficient Asa of the spin axis distance is 

XsA = - cos 6 cos (A + © - a) (3.1-13) 

Combining equations (3. 1-10) to (3. 1-13) gives, together wilh equation 
(2.2-12), the following estimable parameters 

m - [£, 7), Cl 

where 

^ = A - p sin 6 cos ^ ~ x' p sin 6 cos^* - p sin^* cos 6 d6 

t] - d(p cos ^) - P cos^ tan6 d6 - x'p sin^ (3.1-14) 

4 = d(A + © - O') + tan ^ 

The respective coefficients are 

(A) = {1. Asa, ^ 0 _j^) (3.1-15) 

The number of estimable parameters has been reduced to three as compared 
to five (equation (2.2-12)) when 6 is vaiying. 

The estimable parameters for the case of range differencing are 
readily obtained. For a north-south line with stations i and j S 5 mmietric to 
the equator arid for which Pi ~ Pj, the coefficients of (3. 1-15) are identical for 
both stations. The estimable parameters are, therefore, after differencing 
the observation equations: 

pq = [Cj - Cl, -77i, 4j - Cl] 

with 

Cj - Cl = -2 p sin ^ cos 6 d6 

7?j- 7?i = d(p cos^)a - d(p cos^)i ~ 2x p sin ^ (3.1-16) 

4j - 4i = d(Aj - Ai) + 2y tan^ 

is the latitude of the station in the northern hemisphere. Analyzing the 

coefficients (3. 1-15), it is seen that the separation of all three parameters is 
only possible if the observations cover a wide rar^e of the lunar hour angle. 
For observations within a small range of the lunar hour angle one can set 


40 



cos ( A +0 - O') as 1 

sin (A+©-0!)as A+0-a 

ai]d tiae coefficient Asa becomes a constant. Thus, the parameters d6 and x' 
cannot be separated. Because of the geometric observability conditions, the 
useable range of the lunar hour ai^le will be limited. It is, therefore, 
expected that the along-meridian polar motion comi>onent and the-difference 
in spin axis distance (the latter parameter is equivalent to dUj-t if the U-axis 
is located in the meridian of the two stations) will be more affected by ephem- 
eris errors in lunar declination than the across -meridian polar motion 
component. 

The estimable parameters of the east-west line are generally in- 
fluenced by ephemeris errors in declination. Deletii^ terms of the 
order 1/60 and taking 6 constant reduces the coefficients in (3.1-6) 
and (3.1-7) to 

Xg = 2Uj_j sin 6 eos(0-O!)-2Vj_i sin6 sin (©-a) 

S! » =2 cos 6 cos (© - O') 

= 2 cos 6 sin (© - a) 

Combining these three coefficients and taking the relations (2.2-25) 
into account gives the parameters 


m , TT3_i ] 

with 

Cj_i = Vj_idO' + dUj_i+ Uj_ttan6d6 
t?3_i = -Uj_ido' + dVj_| - Vj_itan6 d6. 
The coefficients are 

(A) = (A„ , X ) 

■}~i V}~i 


(3.1-17) 
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Equations (3.1-17) show that the estimates of and will strong y 
depend on ephemeris errors in declination. In actual computation, these 
two parameters correspond to dUj^^ t -and to- dVj-t. " The expressions ’ in 
(3.1-17) can be simplified by introducing a new coordinate system (U), 


U 



whose first axis lays in the "instantaneous" meridian- which goes through the 
center of the points of the east -west line. Equations (3.1-17) then 
become 


Sj-i = Vj-i da + dUj-t 
T)j-t = dVj-i - Vj-itan6 d6 


(3.1-18) 


According to the first equation in (3.1-18), the estimation of 
the earth rotation parameter does not depend significantly on the 
ephemeris uncertainty in declination. The estimates of the earth 
rotation parameter and the across -meridian component of polar 
motion are, therefore, expected to be of the same accuracy. 
However, this equation does not imply that the rotation parameter 
da is obtained independently of crustal motions. The common east- 
west ciustal motions of the stations will be absorbed in da. The 
second equation determines the chord length of the east-west line. 

Its es-timate depends on that of the declination. 


3.2 Station Geometry 


From the previous analysis of estimable parameters, we see 
that each line, consisting of two co-observing stations, yields only two 
accurately determined station parameters. In order to determine all 
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three- orientation parameters at least two lines are needed. Inspect- 
ing the coefficients of equations (3.1-9), it is clear that two stations 
located along the meridian determine polar motion completely, assum- 
ii^ that there is no crustal motion between them. If they are, in 
addition, located symmetrically with respect to the equator, the 
polar motion estimates are virtually independent of an. error in the 
earth rotation parameter da. If the two stations are, however, lo- 
cated on the same parallel, the range difference determines the rota- 
tion parameter da and not polar motion. Lines of north-south and 
east-west directions give the most favorable station geometry since 
the normal matrix of the second adjustment is diagonal for such cases, 
and thus minimizes the estimated correlations between the parameters. 
In view of crustal motion, a clear definition of the terrestrial coor- 
dinate system requires that a third line be added if the other stations 
are located along meridians and parallels. One north-south and one 
east-west line give no over determination for polar motion. 

Therefore, a second north-south line is needed to make polar motion 
free from the effect of individual station motions. Of course, the 
common station motions are still inseparable from the orientation 
parameters. 

For an ideal station distribution, approximate formulas can be 
given for the variance propagation. The configuration consisting of 
two north-south lines of equal lengths, symmetric with respect to the 
equator, and separated in longitude by 90 degrees, is considered 
first. For simplicity, a diagonal variance-covariance matrix of obser- 
vations is assumed in the second adjustment. Under those conditions, 
the normal matrix in equation (2.3-3) is (compare also equations (3.1- 
9)): 
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wf-1 


■FS;-F^ =1 


w!-k 




c?u" 


J-t 


wf-k 

(?v" 

l-lc 


(3.2-1) 


The parameters are the polar motion coordinates x and y. The station 
pairs i,j and k, 1 are eo-observing respectively. Because of the sym- 
metry with respect to the equator, the coordinate differences of the 
third coordinates are equal. 


Wj-t = Wi-k 

Since the two lines are 90 degrees apart, there is a symmetry in the 
coefficient and Ay in (3.1-7) for the respective parameters of 

i-‘ j-i 

the two lines, which, results in 


1 




0?u" 

j-i 




1 



Substituting the two equations above in (3.2-1) and taking the square 
root of the inverse elements gives the standard deviation 


ffu"' cr/ 


Gx S Gy- S' 




l-k 


W 


J-i 


V 




j-t 


i-fc 


(3.2-2) 


If the U-axis is chosen in the meridian of the north-soulh station pair 
i, j, then the variances in equation (3.2-2) are those of the along-meri- 
dian and aeross-meridian component of the estimable parameters. But 
these two variances are characteristically different because the esti- 
mate of the along-meridian component is strongly effected by ephemer- 
is errors in declination. One may, therefore, neglect across -meri- 
dian variance in the denominator of (3.2-2). In doing so, and 
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evaluating the factor pAV^-i, the expressions for the error propagation 
become 


CTt// 

05,/ - CTy/ S 

sin ® 


(3.2-3) 


where $ is the station latitude. Consequently, two north-south lines 
which are separated in longitude by 90 degrees yield the .most accurate 
estimates of polar motion. The accuracy is. the same in all directions 
(circular error ellipse). 

For a single north-south line, symmetric with, respect to the equa- 
tor, the variances are largest and smallest for the along and the across- 
meridian component, respectively. With the U-axis again located in the 
meridian, their accuracies are 


CTy// ^ (3.2-4) 

^ sin^ 


Similarly, an approximate relation for the variance of the rotation 
parameter can be obtained from equation (3.1-18) (east-west line) : 



CT y> 
t j-i 


■cos^> sinu^lL—Ai 


(3.2-5) 


The latter relation holds because the estimate of is insensitive 

to errors in lunar declination. 

It is emphasized that in the subsequent numerical experiment the 
variance of dl/j-'i and dVj-i (equations 3,1-9) will be givdn instead of Cj-i and 
(equations 3. 1-18) since the numerical results are not transformed to the 
special coordinate system (U) and the declination is varied. 
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All these approximate expressions depend on the variance of the 
estimable parameters [X3] of the first adjustment , which, in turn, are 
a function of the number of observations and their distribution over the 
hour angle of the moon. These dependencies are the subject of sub- 
sequent numerical calculation. 

The chords and angles between stations can be computed from the 
estimable coordinate difference in the (U^Vsystem or from the trans- 
formed coordinates in the (U)-system. Both procedures give the same 
results since the orthogonal transformation leaves chords- and angles in- 
variant. Generally, the third coordinate (W-parameter) will set the 
limit to the achievable accuracy. As special cases, the azimuth of a 
north-south line and the chord length of an east-west line are determin- 
able with high accuracy. 


3. 3 Numerical Analysis 


3.3.1 Incrementing the Normal Matrix 


The analysis consists of incrementing the normal matrix based on 
hypothetical observations. The normal matrix of (3.1-2) 


Nx 


■ a;^ PAi 
-AjPAi 


a; PA3- 
A3 P A3 _ 


(3.3-1) 


is formed from the coefficients (3.1-6) to (3.1-8). The respective 
parameters are 

[X] = [Xi i X3I - [dA,d6 i dU/li, dvA» dW/-i, dv'Uu dv/^i, dW^i . . .] 

whereby each new line consisting of two stations adds six parameters 
to the set [ X3]. The inverse of the nonnal matrix is the variance - 
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covariance matrix of the adjusted parameters 

Qx = N-^ 

Particularly, the variance-covariance matrix of the subset [X 3 ] is 

Qx 3 = (A 3 PA 3 - Aa'PAi (AIpAi)'" AIpA3)'^ (3.3-2), 

The a priori and a posteriori variance of unit weight is set equal to unity , 
P is the weight matrix of the range difference observables 

1 

p = iTif I 

corresponding to a single range accuracy of 3 cm. I is the unit matrix. 
Next, the normal matrix (2, 3-3) of the second adjustment 

: (3.3-4) 

is set up and inverted. The orientation parameters [Yg] are 

[YJ = (X, y, da ) 

Finally, the variances of the residuals (2.3-2) are computed by 

Qv = Qx - Qv F (3.3-5) 

33 ® 

3.3.2 Design Characteristics 

Each of the experiments is distinguished by the followir^ design 
characteristics : 

a) number of lines 

b) orientation of lines 

c) length of the lines 

d) observation schedule 

In Table 3.1, column 2, the station configuration for the various 
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Table 3 . 1 

Design Characteristics 


Figure 

Llnes/ObservatlonScheduIe* 

Station Position 

Modification 

3.5 

NS - EW “ NS/A 
165 - 43 - 21) 

NS:A 6 -b-A 3 -i = 90’ EW!A3=80“, A 4 = 20 " 

$3,b=4>^i = 80“ = 


3.6 

as above 

as above 

observation density 

3.7 

NS - EW - NS/B 
(65 - 43 - 21) 

as above 


3.10 

as above 

as above 

50°SAB-e’«Aa_iS90' 

3.11 

NS - EW/A 
(43 - 21) 

NS: 4 > 3-*4 = 80° EW:Ai = 80°, Aa= 20 ” 
4>i = 4>fl = 0 


3.12 

NS - EW/B 
(43 - 21) 

as above 


3.15 

as above 

as above 

(0 . 1 m)® s 1 *^ <!• 4 m)® 

3,16 

as above 

as above 

80°S:4'a-4>4 5:40° 

3.17 

as above 

as above 

azimuth of NS 
0 s: azimuth s: 42° 

3.18 

NS - EW - NS/B 
(65 - 43 - 21) 

NS!A 3 -A 4 = 90'’ EW:Aa=80°, A4=20° 

^>e-s=«>a.i = 80° <E>a = $4 = 0 

20 ° s A 4 S 60° 

3.19 

as above 

Table 3. 2 

*3 and ^4 

3.20 

NS - EW/B 
(42 - 21) 

NS:'1>3-^4=80‘' BWiAi = 80°, Aa=20° 

= 0 

azimuth of EW 
90° & azimuth ^ 48° 

3.21 

conceivable real station distribution (Table 3.3) 



* NS: Stations are located on the meridian. 
EW; Stations are located on the parallel. 




experiments and observation schedules is given. The letters NS and EW 
denote a nortii-south and east-west line, respectively. The first symbol in 
this notation denotes the most western line. It is always a north-south line 
and defines the zero longitude . An exception to this convention is the case of 
Figure 3. 21 where the longitude is zero at Greenwich. Figure 3. 1 shows the 
system of station numbering for the three- and two-line desi^. Thus, the 
lines 56 or 34 define zero longitude, respectively. 


4 

3 

1 

3 


2 

4 


Figure 3. 1 Station Numbering for Three- Line (NS-EW-NS) 
and Two-Line (NS-EW) Design 

Two observation schedules, A and B, were used. In each schedule the 
observation spacing is ten minutes. This time limit was selected for practical 
reasons since it allows one to compute normal points of a certain epoch in 
case of near simultaneous observations . In schedule A, the observations are 
equally spaced throughout the lunar hour angle. The only additional limitation 
is that the altitude of the moon has to be larger than 20 degrees in order to 
avoid disturbing effects of the atmosphere . Schedule B used three hours of 
observations, whereby after the first and third hour there is an interruption of 
one hour. The observation time is arranged symmetrically with respect to 
the lunar transit, i.e., the observations start 2.5 hours before transit and 
finish 2.5 hours after transit. The moon, therefore, has to stay at least five 
hours above an altitude of 20 degrees. 
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0 


Time 


, *-=— f ' 1 1 ■ 

0 1.0 2.0 Time [h] 

Figure 3.2 Observation Schedule 

(The observation, spacing is ten minutes 
within the dotted area.) 


There are two modifications of the station geometry: 

(a) Basic Design : The stations are located exactly along meridians 

(symmetric to equator) and parallels. In case of the three-line 
design, the two north-souih lines are 90 degrees apart. In 
order to demonstrate the influence of the lunar declinational 
position on the variances, the computations are carried out 
for 14 successive intervals, each interval lasting one day. 

In interval 1, the declination is approximately -17.5 degrees, 
between the intervals 5 and 6 it passes zero, and at interval 
13 the maximal declination is reached (Figure 3.3). 

(b) Modified Design : One of the lines changes its position in 

azimuth, length, etc. In all variations the lunar declination 
of interval 1 (-17.5 degrees) is used. 


50 


^max 



Figure 3. 3 Variation in Declination 


3.3.3 A Priori Weighting 

Some consideration has to be given to the a priori weights of the param- 
eters of the first adjustment. Since the coefficients of the parameters duj+i, 
dVj+i, anddWj+i (see coefficients 3. 1-8)- are of order 1/60, these parameters 
are not expected to be of any significance in view of the good approximate, 
coordinates which are available. An initial adjustment was carried out which 
included these parameters. The normal matrix was lU- conditioned. A numer- 
ical inversion with double precision arithmetic on the computer failed. As a 
next step, these parameters were weighted with 

Ou" = Ov'', = CTw" = ± Im 
j+i J+i j+i 

Also the geocentric reflector distance. A, and the third coordinate difference, 
dw^j-i, were weighted constraint with Va = ± 100 m and aw"_ ^ ± 1 m. This 

a priori weightii^ of the station parameters is in accordance with what can be 
expected from present day Doppler positioning. The variances, for the param- 
eter dUj+i, dVj+i, and dWj+i, after inversion of the normal matrix, corres- 
ponded very closely to those in the a priori weights. This, of course, confirms 
that the range difference observables in LLR do not improve these parameters. 
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A final adjustment was made deleting these parameters entirely. The 
resultant variances of the remaining parameters, as wd.1 as the variance of 
unit weight, which was based on .simulated- observations with Gaussian noise 
superimposed, did not differ from the case which included the we^hted 
parameters dUj+i, dVj+i and dW^i, thus making any flirther statistical 
testing unnecessary. 

The parameter dWj-t plays a special role. In all numerical studies 
the a priori wei^t is based on cr ^ = 1 m. The estimated accuracy of the 
declination depends very much on the knowledge of Wj-t. The coefficients 
and .^6 for the north-south line are according to equations (3.1-6) 

Kb “ cos 6 


and 


Aw^_^ “ sin 6 


Since the declination does not change very much during one interval, both 
parameters will be strongly correlated. The numerical effect of the -a priori 
weight of Wj-i on the declination is shown in Figure 3.15. The weighting of 
Wj-i creates another but minor problem. Since Wj-i is a coordinate difference 
in the (U")-system, it is a time varying quantity. The use of the same a priori 
'weight and the same approximate coordinates, e.g., those of the initial epoch 
To, will introduce distortions to the adjustment (too optimistic wei^t). As 
polar motion increases, die once adapted approximate values become increas- 
ingly worse. This problem can be avoided by introducing the adjusted coordi- 
nates of the first adjustment of the previous interval as approximate coordi- 
nates in the subsequent interval. 
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3,3.4 Graphical Rep res entations 

All figures represent the square root of the diagonal element of the 
inverse of the normal matrix (standard deviations) . According to expression 
(3. 3-3), the accuracy for one range difference is 3\/^cm. For other, accu- 
racies ihe graphs change proportionately. For Ihe purpose of checkii^ the 
programming, artificial observations were simulated with a Gaussian noise 
superimposed. The variance of unit weight a posteriori thus obtained fluc- 
tuated around one. But this value is of no real use since it only reflects the 
’’randomness" of the random number generator on the con^uter. It is under- 
stood that with real data, the a posteriori variance of unit weight of the first 

I 

adjustment becomes the a priori variance of unit weight of the second 
adjustment. 

The accuracies of the polar motion coordinates as well as the station 
coordinates are given in centimeters, whereas the a a andCT6 are given in one 
one-lhousandth of an arcsec. All stations are located on a sphere with radius 

R = 6370 km 

The scale of the plots varies among the figures. For each experiment, 
the scale of (ctx, c^y), {au , ay , Uu" . ctv'' 1, and (CTw, o’w" } are, 
respectively, the same. In each case, there is a strong resemblance between 
the standard deviations of the coordinate differences in the (0) -system and 
the (U'^j-system. This is to be expected because of the small degree of 
freedom in the second adjustment. Figure 3.4 reviews the geometric meaning 
of various parameters and related coordinate systems . 
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Figure 3.4 The Coordinate Systems (U) and (U'^) 

The common origin of the coordinate system (U) and (U") is at the 
instantaneous center of mass. For the initial interval at To> tiie two 
coordinate systems coincide per definition. If there were no crustal 
motions, the (U)- system would be fixed to the crust. The body-fixed 
direction of the U-axis is determined by the Greenwich apparent sidereal 
time © and the reflector right ascension a which were used as approx- 
imate (adopted) values in the calculations during the interval To. Note 
that there is no second adjustment for the initial interval To. The 
system (U") is not body-fixed. The third coordinate axis w" coincides 
with the direction of the north celestial pole, and the position of the ^ 
first axis U"is a function of the © and a which are used as approximate 
values in the calculations during a particular interval Ti . In case of , 
crustal motions, the (U)-system is not body-fixed any more but rotates 
slowly by the amount of the common crustal motion component of all 
participating stations.. It is exactly the same motion -which, as mentioned 
earlier, is included in polar motion and earth rotation parameters. 

There is no way of separating this motion from laser ranging at the 
surface only. 

Figure 3.5 shows the basic three -line configuration for observa- 
tion schedule A. Generally, the best estimates are obtained for 6^0, 
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which occurs between intervals 5 and 6 . At that instant, ihe variation in 6 is 
largest; therefore, best parameter separation is possible, and die number of 
observations is largest (longest visibility) . The polar motion accuracy is the 
same in both directions. There is a dependence on the declination, but the 
total variation is very small. The standard deviations of the along-meridian 
components dV^i and dU^s show the ejspected stror^ dependence onithe 
declinations, whereas the across -meridian components dUs-i and dVe-s do 
not show such a dependence . The same dependence is reflected in <y v^_ ^ » 
crug_g andaug_3^, 0 ’Vg_s, respectively. The sudden jump in and 

Cm " results from a change in sign of the correlation coefficient cr i/' 5 as the 
moon passes zero declination. Remembering that the a priori variance of the 
Wj-i parameters is 1 m^, it is seen that only a minor improvement takes 
place. The accuracy estimates of the coordinate differences for the east-west 
linfi CTu ejqiectedly depend very much on the accuracy of Ihe 

declination. .There is, of course, a strong correlation between dV 4-3 and 
dU4-3 since, in the case of only one east- west line, they directly depend on 
the estimated rotation parameter . If the variance -covariance matrix of 
observations for the second adjustment were taken to be diagonal, the corre- 
lation between dV4-3 and dU4-s would be identically 1. Another example of 
the effect of the full variance- covariance matrix of observation is reflected in 
the estimate CTa . Based on ihe corresponding estimable quantities dU 4-3 and 
dV,4-3 , one could again set up an approximate error propagation for a 
diagonal variance- covariance matrix of observation. The geometry is given 
by equations (3.1-9). But such an approximate error estimation results in 
estimates of aS which are too pessimistic since the correlation between the 
observations (estimable parameter of first adjustment) are neglected, .A 
numerical example follows later. 
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*) One unit in abscissa denotes one interval. The lunar declination is varied. 


Figure 3.5 Three-Line Design with Observation Schedule A 
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Figure 3.6 shows the estimated accuracies as a function of the 
number of observations. The basic observation schedule A, which 
included an observation every ten minutes as long as the zenith distance 
of the moon is smaller than 70 degrees, is modified by deleting up to 
80% of the possible observations. Multiplying the number on the horizon- 
tal axis in the graphs by 5 gives the percentage of deleted observations 
from the basic schedule A. The deleted observations are selected randomly 
by use of the random number generator in order to eliminate the 
effect of observation geometry as much as possible. Each computation 
•is made for the same lunar declination (-17.5 degrees). The graphs 
show that a significant deterioration in accuracy occurs after deleting 
50%, and more, of the observations. The accuracy of the across-meridian 
polar motion component (Ty remains unchanged, even for a deletion of 80%. 
The results of this ejqperiment justify the adoption of observation 
schedule B as the standard schedule. 

Figure 3.7 shows the three-line configuration based on observation 
schedule B. There are no new features, detectable. Generally, the accu- 
racy still improves for* zero declination, but not as much as in Figure 
3.5 because of the deletion of observations at the extreme hour ai^le 
in schedule B. 
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*) Multiplying the number in the abscissa by 5 gives the percentage of 
observations which are deleted from schedule A, 


Figure 3.6 Observation Density (Two Lines) 
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*) Multiplying the number in the abscissa by 5 gives the percentage of 
observations which are deleted from schedule A. 

Figure 3,6 (cont’d) Observation Density (Two Lines) 
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*) One unit in abscissa denotes one interval. The lunar declination is varied. 


Figure 3.7 Three-Line Design with Observation Schedule B 
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*)■ One unit in abscissa denotes one interval, lunar declination 
is varied. 

Figure 3,7 (cont'd) Three-Line Design with Observation Schedule B 
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*) One unit in abscissa denotes one interval. The lunar declination 
is varied. 


Figure 3.7 (cont'd) Three- Line Design with Observation Schedule B 
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Figure 3.8 Eigenvalues for the Three-Line Design 


Figure 3.8 shows the eigenvalues for the correlation matrix of the 
observations for the-previous -Ihree-line configuration. Since the TnayiTmim 
and nainimum eigenvalues differ appreciably from unity, one expects the 
correlations of the adjusted estimable parameters [X3] to have a strong 
■influence on the accuracy of the orientation parameters according to the 
inequality (2.3-6). The eigenvalues are slightly affected. by the lunar 
declination'. In Figure 3.9 the accuracy estimates of the parameters of the 
second adjustment are' given for the three-line design using a diagonal 
variance-covariance matrix of observations. -The Cx »- O' y , and Ca have all 
increased. The figures do not include the estimable parameters of the first 
adjustment since they .do not change by this operation. Moreover, the shape 
of the curves has changed significantly. In case of polar motion, one now 
'gets' a. decreased accuracy for zero declination as opposed to an increase as 

* * < ’ i 

in Figure 3. 7. The largest decrease in accuracy occurs for dS . ]h fact, 

CTa is now solely determined from the estimable parameters dU 4-3 and dV^-s , 

t 

which, according to equations (3. 1-17), strongly depend on the ephemeris 
errors in declination. But this apparent dependence betweenCTa and <T& is a 

A 

fallacy. It appears only because the correlations between [X3] were neglected. 
As was shown in equation (3.1-18), the rotation parameter depends only 
insii^ificantly on declination errors .' It Is therefore mandatory to base the 
second adjustment on the full variance-covariance matrix » not only 
because the -variances increase but equally important because the adjusted 
orientation parameters lose their geometric meaning otherwise. 
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*) One unit in abscissa denotes one interval. The lunar 
declination is varied. 

Figure 3.9 (cont’d) Orthogonal Transformation with Uncorrelated 

Observations (Three Lines) 
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In Figure 3.10 the results of var 5 rLng the separation in longitude 
between two north-south lines, are shown. The basic three-line design is 
used. The line 6-5, which defines the zero meridian, remains fixed, 
whereas line 2-1 changes its longitudinal position from 90 degrees to 50 
degrees in steps of five degrees. In ihe figure, only those parameters 
are given which are sensitive to such a change. The y-polar motion 
component, i.e . , mainl y the across -meridian component of the line 6-5 
does not chaise significantly since this line is held fixed. The x-component, 
on the contrary, changes, although not significantly. Therefore, the two 
north-south lines have to be separated only approximately by 90 degrees. 

Figures 3. 11 and 3. 12 show tbe result of the two- line design for 
obseryation schedules A and B, respectively. The only significant difference 
is in the along-meridian polar motion, vhere schedule B gives worse 
accuracy for 6 s 0 because of the neglect of observations at the extreme 
hour angle. The ffy of the across-meridian component indicates only an 
insignificant variation. 
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*) The separation in longitude is reduced by five degress in each step. 
The initial separation in longitude is 90 degrees . 


Figure 3.10 Separation in Longitude for Two North-South Lines 
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*) One unit in abscissa denotes one interval . Th.e lunar declination 
is varied. 


Figure 3.11 (cont’d) Tvvo-Line Design with Observation Schedule A 


71 




























o 

« 


O 

Q — 1 


.(cm) 

10 tl 

J 

1 

no significant 

?• 

O 

o- 

no significant 

n S’" 

variation 

« 2“ 
1 

variation 

,o 

{ L 1 i i 1 1 1 1 1 1 1 

4 

30 

1 . -■■■ 

1 i 1 r i 1 f ! f 1 1 1 


e *1 6 8 10 t? : 

in “ 

2 « 6 8 lO l? ! 


*) One unit in abscissa denotes one interval , The lunar declination 
is varied. 

Figure 3. 12 (cont'd) Two-Line Design with Observation Schedule B - 
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Figure 3, 13 Eigenvalues for Two-Line Design 

In Figure,3.13 the eigenvalues of the correlation matrix of [X 3 ] are 
shown for the two-line design. The maximum eigenvalue is increasing as the 
moon moves from Tninim um to maximum declinations . Figure 3. 14 displays 
the two-line design if the correlation between the observations is neglected. 
All conclusions which were drawn for the tbree-line case, are also valid here. 

The effect of weighting the third coordinate difference on 

the obtainable accuracy in declination is demonstrated in Figure 3.15. 

The following a priori weights are used: 

. 1 ^ 1 

(.Im)® 0 -®// (1.4 mf 

-1 . 

CT6 is nearly directly proportional to the a priori a /, . Polar motion 

. i-i ■ ' 

and rotation accuracies are essentially independent of the weighting. 

• This is important since it confirms that the orientation parameters can 
be obtained accurately even if the third station coordinate is known to 
only ±lm {accuracy of Doppler positioning). Figure 3.16 shows the 
effect of variation in length for a north-south line symmetric with 
respect to the equator! The five experiments are' based on the following 
latitude separation of the stations 

80° s $3 - ^ 2^ 40° 
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*) One unit in abscissa denotes one interval. The lunar declination is varied. 

Figure 3.14 Orthogonal Transformation with Uncorrelated 
Observations (Two- Line Design) 
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Figure 3.15 A priori Weighting of W j - 1 
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Figure 3.15 (cont*d) A 
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The length is decreased in steps of 10 degrees. A long north-south 
line gives obviously a better accuracy in declination and in polar motion. 
According to the approximate formula (3.2-4), and Oyt is inversely 
proportional to the sine of the station latitude. As a whole, the 
deterioriation in polar motion .accuracy as a function of t he latitude 
separation is not very critical, Ar^r separation between 80 and 60 de- 
grees in latitude is acceptable. 

The effect of the variation in azimuth of the north— south line is 
shown in Figure 3.17. The azimuth is changed from 0 degrees to 42 
degrees, in steps of 3 degrees. It is clear from the figures that 
the requirement for an exact north-south line is not very strong. Any 
line, rimning approximately north-south (±15 degrees) will be sufficiento 
In the next three .experiments, the variations of the east-west 
line regarding length, latitude and orientation are investigated. Figure 

3.18 shows the results of an equatorial line with lengths 

60° ^ As - A4 ^ 20° 

The variation occurs in steps of 10 degrees. Only the graphs for 
the east- west line and the rotation parameter are given. Station 3 is 
held fixed at As = 80° and station 4 moves toward station 3. The 
accuracy of the rotation parameter decreases as expected from the ap- 
proximate expression (3.2-5), Generally, a variation in length 
between the limits 60 degrees and 40 degrees is tolerable. In Figure 

3.19 the latitude of the east-west line is varied. Because of the geometri- 
cal constraint imposed by schedule B, i.e. , useful visibility of at 

least 5 hours daily throughout the month, the length of the line has to 
be decreased as the latitude increases. The following positions are 
used: 
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Table 3'. 2- 


Station Variation for East-West Line 
(Variation in Latitude) 


Experiment 


Ai 

Aa 

1 

0 

— 

80 

20 

2 

10 

80 

30 

3 

20 

80 

40 

4 

30 

80 

50 

5 

40 

80 
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•The variation or a agrees with that predicted by the approximate formula 
(3. 2-5). A latitudinal position up to 20 degrees (-20 degrees) for the 
east-west line appears acceptable. 

Finally, Figure 3.20 shows the variation in azimuth of the east- west 
line. The azimuth is varied from 90 degrees to 48 degrees in steps of three 
degrees. Any azimuth up to approximately 20 degrees has an insignificant 
effect on the rotational parameter. 

■' - Figure 3.21 displays the accuracies obtainable for a conceivable real 
station distribution as a fiinction of the lunar declination. 

Table 3.3 

Conceivable Real Station Distribution 


Station 

Latitude 

Longitude 

(1) Texas (McDonald) 

30“ 

256° 

(2) Hawaii 

20“ 

204“ 

(3) Japan 

35“ 

138° 

(4) Australia 

-35“ 

149“ 

(5) Southern Europe 

38° 

15° 

(6) Southern Africa 

-34° 

20° 
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The first four stations are expected to be in operation soon. In this 

( 

ejq^eriment, the zero longitude is at Greenwich. The only new feature in these 
figures is that a a remains small for hi^ lunar declination because'the 
"effective east-west line," i.e., McDonald-Hawaii, is in the northern 
hemisphere’, 
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*) The azimuth varies in steps of three degrees. The initial azimuth is 90°. 
Figure 3.20 (cont’d) Variation in Azimuth of East-West Line 
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*) One unit in abscissa denotes one interval. The Ixmar declination is varied. 


Figure 3.21 (cont'd) Experiment with Realistic Station Distribution 


89 







r* 






*) One unit in abscissa denotes one interval. The lunar declination is varied. 
Figure 3.21 (cont'd) Experiment with Realistic Station Distribution 
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4. SUMMARY AND RECOMMENDATIONS 


In this study the least squares formulation is given for range 
difference observables. The necessary relativistic considerations as 
they relate to high precision lunar laser ranging are given. It is 
pointed out that the earth orientation with respect to the lunar position 
can be parametrized in two equivalent ways. The first method is 
based on the estimable parameters, which are obtained by forming pa- 
rameter combinations (re-parametrization) such that the remaining de- 
sign matrix is non-singular. They are shown to be the coordinates in 
a system (U^^) whose third axis coincides with the celestial pole and 
whose origin is at the instantaneous center of mass. The first axis of 
this system is not strictly fixed to the crust but depends on errors of 
the station clock and the lunar ephemeris. No separate deter- 
mination of the corrections to longitude, time, and right ascension of 
the reflector is possible. The second method uses the estimable quanti- 
ties mentioned above and performs an orthogonal transformation (over 
determined case) so as to result in actual polar motion coordinates and 
in an earth rotation parameter which relate the system (U'') and (U). 

It is necessary to agree upon a standard epoch for which these two 
systems coincide. Although these orientation parameters are sometimes 
referred to as unestimable, they will in no way be inferior to the 
estimable quantities when a standard epoch is fixed. Still, the earth 
rotation parameter is linearly dependent on the correction to lunar 
right ascension. The parameters of the station coordinates appearing 
in the range difference observation equation have been transformed into 
their differences and sums. Because of the earth-moon geometry, the 
coefficients of the coordinate sums are characteristically of the order 
1/60. They were found to be insignificant on the basis of the available 
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approximate station coordinates. The range difference observable will, 
therefore, only determine the coordinate differences, which, however, 
determine" the orientation parameters compretely. Strictly' simultaneous- 
observations are the best since they reduce the significance of ephemeris 
errors. Processing near simultaneous observations requires a good 
ephemeris in order to compute the change in lunar position between 
the twO' epochs of reflection. It is possible to create artificially si- 
multaneous observations using powerful interpolation tools. But such 
investigations should be carried out on real data.- ' . 

An analysis of variance has been based on two types of approximate 
models. The numerical computations have been carried out with a 
model in which the earth rotation during the travel time of the pulse 
was neglected. The accuracy for the range difference was assumed to 
be cm. A second type of approximate model has been set up by 

neglecting terms of the order 1/60 and taking the declination constant 
for the time of one interval. To the approximation of this model' it 
has been shown that the across-meridian polar motion component and 
the earth rotation parameter are independent of errors in declination 
whereas the along-meridian component strongly depends on such errors. 
The analysis was based on ideal station distributions in the form of • 
north-south and east-west lines. Such a design reduces the correlations 
between the parameters. As a basic observation schedule, one observa- 
tion every ten minutes to one and the same reflector was assumed for 
a period of three hours per day. The first and, last observations were 
placed 2.5 hours before and after lunar transit. After one full hour of 
observations, an interruption of one hour was assumed so that the total 
observation span was five hours a day. This condition and the require- 
ment that observations are only made at zenith distances smaller than 

70 degrees put a limit on the station separation. Common visibility 

( 

during the whole month is possible for the following station separations: 
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north-south line (symmetric to equator): 

east-west line: ^ = 0° 

^ = 10 ° 

= 20 ° 

$ = 30° 


< 80° 

A ^ < 60° 
A A- < 60° 
A A < 50° 
A A < 30° 


Each line gives only two estimable parameters accurate enough to 
be useful in deter m i ning the orientation parameters. Usually, the 
accuracy of the third coordinate is modestly accurate at best, because 
of the small change in. declination during one day. Therefore, at least 
two lines are necessary. A three- line configuration is preferable be- 
cause it allows the elimination of individual station motions due to crus- 
tal motion, although the common crustal motion component is absorbed 
by the orientation parameters. Eor a design which includes two north- 
south lines and one equatorial east-west line, the orientation parameters, 
i.e., polar motion and earth rotation variations, can be obtained at 
least with the measurement accuracy. The numerical analysis showed 
that the requirement for lines to run exactly north-south or east-west 
is not very stringent. It is quite sufficient if these two principal 
directions are approximately realized, say, within 10 or 15 degrees. 
Consequently, there is a large degree of freedom for the practical 
realization of such a network. Besides the station geometiy, the 
weather conditions are very important indeed. The final selection of 
station sites should give due considerations to the local climate. Some 
of the stations, which are presently available, already fulfill the geo- 
metric requirements. The stations in Australia and Japan are located 
very ideally to form a north-south line, whereas the McDonald Observa- 
tory in Texas and the station in Hawaii can form the east-west line. It 
is suggested that the method of range -differencing be tested as soon as 
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all of these stations become operational. The missing north-south line can 
be established by mobile laser stations . 

■Gnce accurate orientation parameters become available, they can be 
used as known parameters in long-term single station solutions which will 
provide geocentric coordinates and significant improvements in the lunar 
ephemeris. At that time, it will definitely be possible to reach an accuracy 
level in the ephemeris which allows the analysis of earth core motions 
[Leick, 1978], Yet, an inseparability between the nutations and ephemeris 
corrections exists {Appendix B). This simply demonstrates that the orienta- 
tion of the celestial pole in' space can only be given relative to ’the motion of 
the moon in case of lunar laser ranging. 
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APPENDIX A 


Generalized Inverse Solutions and Estimability 


The two step procedure to solve for polar motion and earth rota- 
tion described in the previous sections can also be formulated as a 
one step procedure. Such a re-formulation cannot effect the outcome 
of the adjustment, but it helps in demonstrating commonly used termin- 
ology. 

Consider the equation (2.2-30) and the following equation which 
gives the general model for range differences as 


BV + (Ai Ag As) 


■Xi“ 

Ys 

-Y3- 


+ W = 0 


and which includes all the parameters explicitly. This adjustment 
model, in which each equation contains two range observations together 
with the parameters, had to be used because of the finite velocity of 
the light. For the purpose of this appendix, we can limit ourselves to 
the simplified model expressed by equations (3.1-6) to (3.1-8) where 
the earth rotation during the travel time of the pulse is neglected. 
These simplifications led to the model with observation equations 


■X, 


V ~ (Ai Ag As ) 


+ L 


LYs-I 


(A.l) 


where 



(Sy/ 

•^3 ” (-^U > -^V > J ^V J -^W ' 

3-i 3-1 3-1 3t-i «-l 

V are the residuals and L are the observations. 

The equation (A. 1) is written as 


V = AY + L 


(A. 2) 


with 

A = (Ai Ag A 3 ), (A. 3) 

The design matrix A has a rank deficiency of three. -The weight 
matrix of the observations L is 


P = — I (A. 4) 

3 

. CTo 

The most general solution to this over determined system is the unique 
minimum norm least-squares solution [Rao and Mltra, 1971, p. 51] 
which satisfied the conditions 


and , 


The solution is 


V’' pv = minimum 
QY = minimum 


Y = ApcjL. 


(A. 5) 
(A. 6). • 

(A. 7) 


A pQ is called the minimum Q-norm P-least squares inverse. It fulfills 
the following four conditions fiao and Mitra, 1971, p. 52] if P and Q 
are positive definite (p.d. ) matrices: 


AApQ A = A 


+ + 

A p q A A p q 


= A 


.+ 


PQ 


(AA pq) P - PAA pq 
(A pq A)^ Q == Q A pq A 


(A. 8) 
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If Q is positive semi-definite (p.s.d.), the first two conditions in (A. 8) 
axe replaced by ' 


PAA^ppA = PA 

+ + + 
QAp(jAA pQ =QApQ 


(A. 9) 


The minimum norm least squares inverse A p q is given by the- G-matrix 
of equation (2.3-5) if we replace the matrix M'^' by P of (A. 4) in the 
expressions (2.2-27), (2.3-3) and (2.3-4). Denoting the thus obtained 
matrix by gives 



(A. 10) 


with 



O 


O O 

O O 


-X- 



(A. 11) 


That the generalized inverse indeed fulfills all conditions in (A. 8) and 
(A. 9) can be verified by straightforw^d matrix multiplication. It is, 
therefore,- formally established that the two step solution is identical 
to a minimum norm least squares solution. 

In the main bo(fy of this study, we called the parameters [X^Xs], 
which are the ephemeris corrections and the coordinates in the (U^'")- 
system, estimable parameters. This was done in order to underline 
that the corresponding design matrix was of full rank. According to 
Eao [1965, p. 224], ^ linear parametric fimctions are estimable, if 
and only if the rank of the design matrix is full, i. e. , if the rank is 
equal to the number of parameters to be solved. It is, therefore, 
always possible to find estimable parameters simply by inspecting the 
coefficients of the design matrix for linear dependencies and combining 
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the corresponding parameters to form new but estimable parameter 
combinations. If one does not change the design matrix, i.e. , leaves 
thie linear dependent columns included, there is the following condition 
for unbiased estimation [Rao and Mitra, 1971, p. 139]: A parametric 
function p^Y is unbiasedly estimable by a linear function of L (the 
observations) if and only if 

p^ A A = p^' (A. 12) 

or 

p''‘(A^'A) A^A = p*'^ (A, 13) 

where A is any generalized inverse which fulfills 

A a” A = A (A. 14) 

From (A. 12) and (A. 14) it is clear that for 

= A A 

the linear parametric function p’^Y is unbiasedly estimable. It was 
shown that A^pq fulfills the condition (A. 14) which is a special case of 
the first condition in (A. 8). With equations (A. 3) and (A. 10), one 
obtains the estimable . parametric function as 

p^Y = a\qAY = 

Using the expression (2.3-2) for [Ya], i.e., 

Ya = -F" Ys + Xs 

the equations (A. 15) can be* rewritten as 


I O 
O I 

O O 


(FS\ FT)-^FSa 

Xa , Xg, 

I - F^(FS^ F^)'^FS- 

3 



“xf 


Ys 

d 

_’^3_ 


(A. 15) 
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I o 

p^Y = O (FS| fV^ fs| 

O I - F^ (F F^)‘^ F S| 

These equations express the parametric function. p^Y in terms- of the 
estimable parameters [Xi, X3] . Corc^aring (A. 16) witii the least 
squares estimates in equation (2. 3-5), the parametric function becomes 

Xi 

p'y = Y2 

A 

.Ys 

It is recognized that both procedures, i. e. , finding the estimable linear 
parameter combinations and performing a second adjustment, or using 
the formalism of generalized inverses, leads to the same result. 

Finally, the present solution is compared with what is sometimes 
referred to as inner constraint or "pseudo inverse solution." At the 
outset, it is underlined that so fer only the subset [Y 3 ] takes part in 
the minimization of the second adjustment. The parameters [X^], 
i.e., the lunar declination and the geocentric reflector distance, are 
entirely independent of the definition of the coordinate system (U). The 
declination refers to the (U^^) -system. Therefore the parameters [Xx] 
are not a subject of any constraint whatsoever in determining polar 
motion. 

Pope [1971] discusses the use of the Null space in solving singu- 
lar geodetic- systems. It is understood that polar motion is assumed to 
be known in that context. The singularity results from a lack of defini- 
tion in shift, possibly in scale, and in a single rotation. The 
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singularities are eliminated by incorporating a similarity transformation 
on the approximate station coordinates in the form of a constTaint. He 
points out that any basis E of the design matrix is suitable in order 
to obtain an inner constraint solution. Adding the constraint 

E’’Y = O (A. 17) 

where E fulfills the condition 

AE = O (A. 18) 

to the normal equations gives an alimented non- singular normal matrix 
whose inverse is 

N 
E^ 

N is the pseudo inverse. It fulfills the conditions (A. 8) with 
P = I and Q = I. The parameters are 

Y = A^ PL = (A^ PA)"^A^ PL . 

. ^ A p ^ L » 

The latter equality is readily proven by use of the properties (A. 8). 
Since the norm matrix, Q = I, is an identity matrix, all 
parameters participate in the minimization, i.e., 

Y"^ Y = Xl Xi + Y ^3 Yg + Y^ Ya = minimum 

The matrix E is readily derived from our previous work. This is most 
easily seen by looking at the relations involved in the rank factorization 
theorem which was used repeatedly in section 2. Graybill [1961, sec- 
tion 11.2.3] , in what he refers to as re-parametrization, gives the 
following relations. Consider the model of observation equations as 
in (A. 2) 

V = AY + L 


i 


E 

.0 


N 


'X ttiT 


(E" 


E(E’' E)'^ 
O 


(A. 19) 
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with an m X m normal matrix 

= (A’’ A) 

of less than fiill rank 

R(N) = m - s = u 

The weight matrix of observations is assumed to be equal to the identi- 
ty matrix for the present purpose, which can always be accomplished 
by transformation. According to well-known theorem in linear algebra, 
there exists non-singular matrix 


“ (b^u bE«) 

such that 


M^(A^ M = 


’uS’’ (A^ A)S„ 
O 


O 

O 


(A. 20) 


where the non-zero submatrix in (A . 20) is of size and rank u. The 
relation (A. 20) implies 


E^{A'^A)e = O 


which in turn implies 

A E = O 


(A. 21) 


Equation (A. 20) also implies that R(AS) = u, i. e. , the product is of 
full rank. E spans the Null space of the design matrix. The original 
observation equation (A. 2) can be rewritten as 

V = AY + L = AMM'^Y + L 


Partitioning the inverse of M by 

= 


_H. 


leads to 
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V = ASHY + AEHY + L 


- ASHY + L ■ ' ■ (A. 22) 

The last equality follows from (A-21). Denoting the product of A and 
S by D, 

D = AS (A. 23) 


one obtains the rank fectorization theorem by comparing (A. 22) and 
(A. 2) as ' - - 

A = DH (A. 24) 


where D has full rank. The estimable parameters result from the 
non-estimable parameters by the transformation 

X = HY (A. 25) 

According to (A, 22), the matrix E is now readily obtained. The design 
matrix (D) of the estimable parameters [X] was given in Section 2.2.2 by 

D = AS = (A 1 A 2 A 3 ) S = (4 A 3 ) (A. 26) 


Equation (A. 26) determines the matrix S as 


S = 


I 

O 

O 


O 

0 

1 


Since 

MM’\= SH + EH = I 
substitution of (2.2-31) for the matrix H gives 


" I 

0 

o' 


'0 

0 

o' 

0 

0 

0 

- 

0 

I 

0 

0 

FT 

I 


0 

-F^ 

0 _ 


which leads to the following identities; 
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r O n 


and 



(A. 27) 


H 5 (o I O) 


The matrix E, which spans the basis of the Null space of. the design 
matrix, can be easily found directly by solving the ^uation (A. 21). 

The identity (A. 27) can be verified using the coefficients (3.1-6) to 
(3.1-8). 

It is concluded that the procedure of finding the basis of the de- 
sign matrix and applying it straightforwardly as a constraint EY = 0 
leads to an equally weighted minimization of the squares of all param- 
eters. This is an undesirable procedure for the case of range differ- 
ences to the moon. First, the [Xi] parameters (geocentric reflector 
distance and declination) should not be included in the minimization. 
Second, the parameters [Ya], i.e. , the station coordinate differences 
and sums, should not be minimized based on equal weights, because • 
some of those parameters are weakly determined. Third, the inclu- 
sion of the orientation parameters, in particular polar motion, in the 
minimization adds a time varying component to the station coordinates. 
The resultant station coordinates are neither crustal fixed (even if no 
crustal motions occur) nor are they components in the celestial system. 
Since the primary concern of this study is not only the determination 
of chords and angles, but also of finding unique orientation parameters, 
the pseudo inverse solution is not pursued any further. 
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APPENDIX B 


Ephemeris Errors Vs. Corrections in Nutation 

Polar motion modelling refers to the celestial pole (C). It was 
modelled as a constant per day, i.e. , the mean of the progressive 
Chandler motion per day. The celestial pole, per definition, has no 
periodic diurnal motion relative to the crust. The adjusted parameters 
in declination and right ascension (the latter is linearly dependent on 
correction to time) will give the corrections to the l unar ephemeris in 
the celestial system (X). These corrections contain two lypes of errors 
which are not separable immediately. One error source is that the 
adopted set of nutations, whichever set one uses, is unlikely to describe 
the celestial pole (C) in space accurately, since any nutation set is de- 
rived from a hypothetical earth model. The adopted pole (adopted 
nutations) will have a nearly diurnal periodic body- fixed motion. The 
second type of error denotes actual ephemeris errors. Ephemeris 
errors can have several origins, such as errors due to truncation, 
errors in the constants, or even programming errors. The latter 
error source is particularly suitable to demonstrate the consequences 
of the non-separability of the two errors. Assume that during the 
programming of the ephemeris one term was forgotten. Not knowing 
this, one interprets the adjusted corrections in declination or right 
ascension as an error of the adopted set of nutations . If this set of 
nutations is compared with another experimental set as derived from, 
say, VLBI, there will be a discrepency just equal to the forgotten 
ephemeris term. 
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I'he polar motion which refers to the pole as defined by the adopted 
set of nutations is decomposed as 


X = Xc + x(t) 
y = yc + y(t) 


(B.l) 


Xc and yc are constant for the time of one interval. They are the components 
of the celestial pole which have always been referred to in Section 2. The 
diurnal components are modelled with the same frequencies as have the nuta- 
tions ! It is assumed that the diurnal component can be totally accounted for 
by changing the coefficients of the adopted set of nutations. 

Equation (2.3-13) in [Leick, T978] gives for the diurnal polar motion 
terms the expression 


x(t) + iy(t) 


-ilZ -A-j 
} 


-i(G]V[ST + Attj) 

© 


= [-Aj sin (GMST + Attj) - i Aj cos (GMST + Aoij)] 

J 

where Aj is the coefficient, A^j the nutation argumoitfor the nutation j, and 
GMST stands for Greenwich mean sidereal time . The y-axis is taken positive 
along A = 90°, whereas in the preceding part of this report y was taken positive 
along A '= 270°. Chai^ng the sign of y in the above equation, the residual 
diurnal polar 'motion is modelled as 


x(t) + iy(t) = ^ [-dAj sin (GMST + Aa^) + idAj cos (GMST + Aaj) 

^ ^ (B.2) 

In [Leick, 1978, Section 2.4], expressions for the change in declination and 
right ascension were given as a function of the nutation coefficients and fre- 
quencies. With Oi being the right ascension, the change in declination is 
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SaaopM - 6c = -Re (iSAje 

) 

= S -Aj sin (AQ!j + a) 
i 

Thus, the change in declination due to an error in the nutation coefficients 
is 

d6(t) = S -dAj sin(AQ!j + a) (B.3) 

The change in right ascension is given by the same reference 

°aaopted-o^c = Im(iL Aj ) tan a ~ ZAj cos(Ao!+ oj tan 6 

^ s 

Thus 

do!(t) = SdAjCos (AQ!j + oi)tan6 (B-4) 

Note that the sign in (B.3) and (B.4) is always in the sense "adopted 
minus C." 

The coefficients for polar motion, declination and right ascension 
are given by equations (3.1-3) and (3.1-4). After some lengthy alge- 
braic manipulations, we find the following relations 

A^ x(t) + 'Ky y(t) + Ae d6(t) + Aa da(t) = 0, (B.5) 

with 

A ^ 

a A + ©-a 


The first two terms are identified as those terms in the observation 
equation which relate to diurnal polar motion. They are not included 
in the models discussed in the main body of this study. The linear 
relation (B.5) shows that errors in the nutations, i.e. , the adopted set 
of nutations does not describe the direction of the celestial pole (C),will 
be absorbed in the daily adjusted declination and right ascension param- 
eters. The nutation errors will be inseparable from any errors in lunar 
ephemeris. It is important to realize that the actual polar motion 
coordinates, that is, the constant part (neglecting the progressive 
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Chandler motion) is obtained uniquely. The adjustment station positions 
are independent of the relation in (B.5). If one were to compare the, 
path of the celestial pole C with respect to the crust as obtained 
feom LLR and, say, VLBI, there should be full agreement regardless 
of whether a "programming error" occurred during ephemeris imple- 
mentation, or whether no such error’ occurred. The comparison can 
be made with the respective polar motion coordinates provided both 
parties select the same "zero point" for counting polar motion. 
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